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ABSTRACT 


The  propagation  of  stress  waves  due  to  a  point  type  excitation  In 
the  form  of  a  sinusoidal  pulse  In  an  Infinite  medium  with  Inclusions  hav¬ 
ing  different  properties  Is  studied.  The  solution  Is  carried  out  using 
the  boundary  element  method  In  the  frequency  domain  with  a  Discrete  Four¬ 
ier  transform.  The  Inclusion-medium  interfaces  are  discretized  using  a 
constant  element  which  assumes  a  uniform  stress  and  displacement  field 
over  the  element.  Studies  were  conducted  primarily  with  a 
two-dimensional  plane  strain  model  but  some  were  also  performed  In  the 
three-dimensional  case,  focusing  on  the  attenuation  characteristics  and 
the  velocity  of  the  wave  In  terms  of  the  arrival  time  for  both  the  free 
field  and  the  case  with  Inclusions.  Results  are  presented  In  the  form  of 
a  dimensionless  displacement  and  arrival  times  at  the  target  under  con¬ 
sideration.  With  a  point  excitation,  as  used  in  this  study,  the  free 
field  attenuation  follows  the  geometrical  damping  law  for  both  the  two 
and  the  three-dimensional  cases,  except  at  distances  In  the  neighborhood 
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of  one  wavelength  or  closer,  where  a  more  complex  pattern  of  waves  Is 
developed.  Some  comparative  studies  were  also  performed  using  the 
finite  element  method  by  Integrating  the  equations  of  motion  In  the  time 
domain  with  an  explicit  Integration  scheme.  The  results  were  similar  to 
those  of  a  system  with  a  small  amount  of  damping.  Cases  with  different 
sizes,  shapes  and  properties  of  Inclusions  as  well  as  different  arrange¬ 
ments  of  an  Inclusion  cluster  were  studied  and  their  effects  on  the 
attenuation  and  velocity  of  the  wave  were  investigated.  Results  show 
that  while  the  effect  on  the  amplitude  and  the  arrival  time  depends  on 
the  wavelength  and  the  relative  size  of  the  Inclusions,  the  attenuation 
rate  Is  practically  unchanged.  The  arrival  time  of  the  wave  Is  affected 
more  by  a  single  Inclusion  than  a  cluster  of  Inclusions  with  the  same 
total  area. 
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Uo  free  field  displacement  vector. 

U(w)  frequency  response  for  a  unit  concentrated  harmonic 
excitation. 

Z  argument  for  the  modified  Bessel's  function, 

a  constant 

*  unit  weight,  Euler1  s  constant. 

T  symbol  representing  boundary. 

6,j j  Kronecker's  delta. 

&  Dirac  delta  function. 

e  small  radius  of  a  semi  circle  or  a  hemisphere  In  the 

boundary  element  formulation. 

unit  normal  vector  component. 

X  wavelength. 

v  Poisson's  ratio. 

p  mass  density. 

w  circular  frequency. 

Q  circular  frequency,  symbol  representing  a  domain, 

x  function  deflnded  In  text. 
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function  deflnded  In  text. 

length  of  one  dimensional  element,  least  dimension  of 
two  or  three  dimension  elements,  subscript  Indicating 
point  on  the  boundary  element. 

superscript  Indicating  the  virtual  state. 

a  cap  Indicating  the  transformed  domain  parameter. 
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CHAPTER  1 

INTRODUCTION 


The  study  of  elastic  stress  wave  propagation  has  regained  momen¬ 
tum  In  recent  years  because  of  Increasing  needs  for  accurate  Information 
In  the  field  of  seismology,  earthquake  engineering  and  soil  dynamics. 
Recent  developments  In  electronic  Instrumentation  have  made  the 
detection  of  small  amplitude  motion  due  to  stress  waves  possible,  and  the 
experimental  work  In  this  area  can  now  cope  with  the  advances  that  had 
been  made  on  the  theoretical  side.  This  has  resulted  In  an  Increasing 
number  of  papers  on  both  the  experimental  and  the  theoretical  aspects  of 
the  problem.  Another  reason  for  this  change  Is  the  development  In  the 
field  of  material  science  with  the  application  of  the  stress  wave  tech¬ 
nique  In  the  Investigation  of  material  properties.  Finally,  the  appli¬ 
cation  of  seismic  techniques  in  geophysics  and  the  study  of  the 
propagation  of  stress  pulses  of  large  amplitude,  which  Is  particularly 
Important  from  the  engineering  point  of  view,  make  the  subjects  even  more 
Interesting. 

In  dealing  with  wave  propagation  In  soil  deposits  where  the 
state  of  stress  varies  depending  on  the  geological  conditions,  the 
Information  on  the  effect  of  the  general  state  of  stress  on  the  velocity 
of  waves  Is  very  Important,  but  the  experimental  results  available  on 
this  area  are  still  limited.  A  research  program  was,  therefore,  con- 
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ducted  at  the  Geotechnical  Engineering  Division,  University  of  Texas, 
Austin,  to  Investigate  these  effects  at  low  levels  of  strain.  A  large 
scale  trlaxlal  testing  device,  which  Is  a  reinforced  steel  cube  designed 
to  hold  a  cubic  soil  sample  measuring  7  ft.  on  each  side,  was  used.  The 
waves  are  generated  along  the  principal  axes  through  one  of  three  excita¬ 
tion  ports  located  on  three  walls  of  the  cube.  Accelerometer  packages 
are  embedded  In  the  soil  along  each  principal  axis  of  the  cube  to  monitor 
the  waves. 

The  first  motivation  for  the  present  work  arose  from  the  need  to 
Investigate:  (1)  the  effect  of  the  much  stlffer  (almost  rigid)  acceler¬ 
ometer  packages  mentioned  above  on  the  characteristics  of  the  waves  pro¬ 
pagating  through  the  soil,  and  (2)  how  well  the  motion  of  the 
accelerometer  package  Itself  represents  the  soil  motion.  The  original 
Investigation  was  conducted  using  the  finite  element  method  with  an 
explicit  Integration  scheme,  which  was  appropriate  for  this  particular 
problem  because  of  the  finite  domain.  The  results  of  the  first  Investi¬ 
gation  were  presented  in  Stokoe  et  al  (1980).  Further  studies  were  also 
performed  on  the  same  problem  using  the  same  numerical  model  to  Investi¬ 
gate  the  effect  of  different  boundary  conditions,  different  types  of 
excitation  and  different  wavelengths  (or  frequencies  of  excitation). 


The  results  of  these  studies  were  presented  In  Suddhlprakarn  (1983). 

There  are  Increasing  Interests  In  the  problem  of  diffraction  and 
scattering  of  elastic  waves  by  obstacles  or  Inclusions  in  recent  years, 
because  of  its  wide  range  of  applications  in  engineering,  geophysics, 
acoustic  and  military  applications.  It  was,  therefore,  decided  to 
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Investigate  further  the  more  general  area  of  wave  propagation  In  a  medium 
with  Inclusions  having  different  properties.  Instead  of  a  finite  domain 
representing  a  large  trlaxlal  cube,  as  used  In  the  previous  studies,  a 
full  space  was  considered.  The  main  reason  for  this  was  to  eliminate  the 
boundary  effects  and  to  concentrate  only  on  body  waves. 

In  this  work,  the  boundary  element  method  Is  used  and  the  prob¬ 
lems  are  solved  In  the  frequency  domain  with  the  application  of  the  Four¬ 
ier  transform.  The  numerical  model  consists  of  an  Infinite  space  with  a 
unit  point  excitation,  Inclusions  and  target  points  where  the  results 
are  calculated.  The  target  points  can  also  be  located  at  the  Inclusion 
surface.  The  number  of  Inclusions,  their  shapes  and  their  arrangements 
can  be  varied  arbitrarily  as  long  as  the  number  of  boundary  elements  used 
does  not  exceed  the  maximum  limit.  The  properties  of  a  medium  and  the 
Inclusions  can  also  be  selected  to  handle  a  wide  range  of  situations. 
Both  three-dimensional  and  two-dimensional  plane  strain  models  are 
developed.  The  formulation  based  on  the  boundary  element  method  and  a 
frequency  domain  analysis  Is  presented  In  Chapter  4. 

Some  results  are  also  obtained  using  the  finite  element  model  to 
compare  with  those  from  the  boundary  element  method.  The  finite  element 
model  from  those  of  the  early  Investigations  Is  modified  to  handle  an 
Infinite  space  problem  as  considered  In  this  work.  The  formulation, 
which  was  not  shown  In  the  early  work,  Is  presented  In  Chapter  3. 

In  Chapter  5,  a  preliminary  Investigation  on  the  effect  of  some 
Important  parameters  Involving  both  numerical  models  are  presented.  In 
the  finite  element  model,  the  investigation  Involves  the  size  of  the 


domain  necessary  to  represent  an  Infinite  medium.  In  the  boundary  ele¬ 
ment  model,  most  of  the  Investigation  focuses  on  the  Fourier  transform 
parameters. 

Chapter  6  shows  the  effects  of  Inclusions  of  different  sizes  and 
properties  on  the  waveform  and  the  attenuation  characteristics  of  the 
wave.  In  Chapter  7,  the  arrival  time  of  the  wave  and  the  effect  from 
Inclusions  having  different  sizes,  shapes,  and  properties  are  presented. 
The  effect  of  an  Inclusion  cluster  versus  a  single  Inclusion  are  also 
discussed.  Conclusions  from  the  results  In  Chapter  6  and  Chapter  7  are 
presented  In  Chapter  8. 
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CHAPTER  2 


BACKGROUND 


2.1  WAVE  PROPAGATION 

The  study  of  elastic  wave  propagation  originated  from  the 
attempt  to  explain  the  true  nature  of  light,  which  was  thought,  In  the 
early  nineteenth  century,  to  be  the  propagation  of  a  disturbance  In  an 
elastic  aether.  After  the  research  of  Fresnel  (1788-1827)  and  Young 
(1773-1829)  on  the  polallzatlon  of  a  light  beam,  It  was  concluded  that 
light  could  be  explained  only  by  transverse  waves  In  an  elastic  aether. 
This  conclusion  paved  the  way  for  further  studies  of  elasticity  and  wave 
propagation  by  many  great  mathematicians.  Navler  (1785-1836),  who  was 
the  first  to  Investigate  the  general  equations  of  equilibrium  and 
vibration  of  elastic  solids,  presented  a  molecular  theory  of  an  elastic 
body,  giving  an  equation  of  motion  In  terms  of  the  displacement  of  the 
material  point  (molecule). 

In  1822,  Cauchy  (1789-1857)  developed  what  Is  known  today  as  the 
"Mathematical  Theory  of  Elasticity",  Including  the  stress  and  displace¬ 
ment  equation  of  motion  which  agreed  with  Navler* s  in  some  aspects.  In 
1828,  Poisson  (1781-1840)  succeeded  In  solving  the  differential  equation 
of  motion  for  an  elastic  solid  by  decomposing  the  displacement  Into  an 


irrotatlonal  and  equivolumlnal  part,  with  each  being  a  solution  of  a  wave 
equation.  This  led  to  the  conclusion  that  an  elastic  disturbance  was  in 
general  composed  of  two  types  of  fundamental  waves,  the  longitudinal 
(Irrotatlonal)  and  transverse  (equlvoluminal )  waves. 

Further  Important  Investigations  made  earlier  on  the  propagation 
of  elastic  waves  were  carried  on  by  several  Investigators  such  as  Ostro- 
grasky.  Green,  Stokes,  Lame,  Clebsch,  and  Chrlstoffel.  The  history  of 
the  studies  on  elastic  waves  and  theory  of  elasticity  is  presented  In 
detail  by  Love  (1944). 

In  the  transition  period  between  the  end  of  the  nineteenth  cen¬ 
tury  and  the  beginning  of  the  twentieth  century,  the  study  of  waves  in 
elastic  solids  regained  momentum  due  to  applications  In  the  field  of  geo¬ 
physics.  There  are  several  contributions  in  this  period  especially  In 
the  area  of  wave  propagation  In  an  elastic  half  space.  In  1887,  Rayleigh 
discovered  his  now  well-known  surface  wave,  which  Is  a  very  significant 
type  of  wave  Involving  a  traction  free  surface.  In  1904,  Lamb  was  the 
first  to  study  the  propagation  of  a  pulse  In  an  elastic  half  space.  This 
problem  Is  of  primary  Importance  In  the  field  of  seismology,  and  was  fur¬ 
ther  studied  by  many  researchers  such  as  Nakano,  Lapwood,  Sobolev,  Smir¬ 
nov,  Stoneley,  and  Cagnlard. 

Since  World  War  II,  the  research  on  elastic  wave  propagation  has 
expanded  at  an  Increasing  rate  due  to  the  need  for  Information  on  earth¬ 
quake  phenomena,  nuclear  explosions  and  the  performance  of  structures 
under  high  rates  of  loading.  A  survey  of  most  of  the  contributions  In 
this  field  until  1965  was  cited  In  Mlklowltz  (1966).  Ewing  et  al  (1956) 
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gave  a  comprehensive  treatment  of  elastic  waves  with  extensive  refer¬ 
ences  on  most  topics.  More  recent  comprehensive  books  on  general  elastic 
wave  propagation  have  been  written  by  Achenbach  (1973),  Eringen  and 
Suhubi  (1975),  and  Miklowitz  (1978). 

2.2  DIFFRACTION  OF  ELASTIC  WAVE 

The  term  diffraction  was  given  by  Grimaldi  (1618-1663)  to  the 
phenomenon  of  a  light  beam  bending  slightly  while  passing  the  edge  of  an 
aperture.  Nowaday,  this  term  is  applied  to  a  phenomenon  of  wave  propa¬ 
gation  when  the  wave  rays  deviate  from  rectilinear  paths,  which  cannot  be 
interpreted  as  reflection  or  refraction.  The  definition  of  a  scattered 
wave,  on  the  other  hand,  Is  not  as  clear.  Rayleigh's  definition  Is  that 
the  scattered  wave  Is  tha  difference  of  the  total  wave  field  observed  in 
the  presence  of  an  obstacle  and  the  incident  wave  (or  the  free  field) 
which  means  that  the  scattered  wave  may  contain  the  reflected  wave  or  the 
combination  of  the  refracted  and  diffracted  wave  depending  on  the  posi¬ 
tion  of  the  target.  Outside  the  field  of  molecular  physics,  the  term 
scattering  and  diffraction  are  often  used  to  describe  the  same  wave  phe¬ 
nomenon.  When  the  diffracted  wave  forms  the  major  part  of  the  scattered 
wave,  as  In  the  case  of  an  obstacle  with  sharp  edges,  the  term  "diffrac¬ 
tion"  Is  normally  used.  On  the  other  hand,  when  the  diffracted  part  of 
the  wave  Is  less  significant  as  In  the  case  of  an  obstacle  without  sharp 
edges,  the  term  "scattering"  Is  preferred. 


In  parallel  to  the  development  of  the  elastic  theory  of  light, 
the  studies  on  diffraction  and  scattering  phenomena  was  also  developed. 
The  first  attempt  was  made  by  Stokes,  who  studied  the  diffraction  of 
light  by  an  aperture  In  a  screen,  and  Clebsch,  who  studied  the  reflection 
and  transmission  of  light  through  lenses.  A  few  years  later,  Rayleigh 
discussed  the  scattering  of  light  by  very  small  particles  In  relation  to 
the  wavelength  and  explained  why  the  sky  Is  blue,  Rayleigh  was  also  the 
first  to  study  the  scattering  of  sound  waves  by  a  spherical  obstacle  hav¬ 
ing  different  density  and  compressibility  than  the  surrounding  air.  He 
considered  plane  as  well  as  spherical  waves  and  the  results  confirmed  his 
law  of  scattering  known  as  Rayleigh's  scattering,  which  involved  a  small 
scatterer  and  large  Incident  wavelengths. 

There  were  several  further  studies  of  the  diffraction  and  scat¬ 
tering  on  both  electromagnetic  and  elastic  solid  waves  by  various  Inves¬ 
tigators.  The  Interest  In  the  diffraction  of  waves  In  an  elastic  solid, 
however,  had  declined  for  some  time.  After  the  paper  published  In  1927 
by  Sezawa  of  the  Tokyo  Imperial  University  on  the  scattering  of  an  Inci¬ 
dent  congressional  wave  by  circular  and  elliptical  cylinders  as  well  as 
by  a  sphere,  numerous  publications  followed,  all  of  which  dealt  only 
with  a  sinusoidal  Incident  wave,  l.e.  with  a  steady  state  response.  The 
problems  of  the  transient  response,  on  the  other  hand,  have  been  Investi¬ 
gated  only  In  recent  years  with  the  first  contribution  by  Baron  and  Mat¬ 
thews  (1961)  who  calculated  the  stress  response  at  the  surface  of  a 
cylindrical  cavity  by  a  step  congressional  wave.  Pao  and  Mow  (1973)  pre¬ 
sented  a  comprehensive  review  of  literature  and  treatment  on  the  subject 


focusing  on  the  dynamic  stress  concentration.  Some  solution  methods  and 
results  on  the  diffraction  of  elastic  waves  by  circular  and  spherical 
scatterers  with  emphasis  on  the  transient  problems  are  also  presented  In 
Mlklowltz  (1978). 

Due  to  the  Increasing  complexities  of  the  geometry  of  the  prob¬ 
lem  together  with  the  development  of  powerful  digital  computers  during 
the  last  decade,  more  research  has  been  oriented  toward  the  application 
of  numerical  models  to  solve  the  problem  of  wave  propagation  and  diffrac¬ 
tion. 

2.3  BOUNDARY  ELEMENT  METHOD 

Most  of  the  problems  In  the  field  of  engineering  and  science  are 
usually  formulated  by  considering  the  behaviour  of  an  Infinitesimal  ele¬ 
ment  of  the  system  with  an  assumed  relationships  between  the  variables 
Involved.  This  leads  to  a  set  of  differential  equations  that  describes 
the  nature  of  the  system.  It  had  long  been  realized,  especially  for  the 
case  of  partial  differential  equations*,  that  general  solutions  are  not 
so  useful  as  the  solutions  for  specific  problems  where  the  Initial  and 
boundary  conditions  are  given.  The  solutions  of  such  Initial  value  or 
boundary  value  problems  has  been  of  major  concern  to  Investigators  for 
many  years. 

The  majority  of  the  practical  problems  often  Involve  complicated 
boundary  conditions  and  parameters  which  are  not  tractable  with  analyt¬ 
ical  solutions.  Numerical  methods  become  then  the  only  feasible  means  to 
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obtain  reliable  results.  The  most  widely  used  numerical  methods  at  pres¬ 
ent  manipulate  directly  the  original  differential  equations  using  vari¬ 
ous  approximation  techniques.  These  approximations  can  be  done  either 
by  replacing  the  differential  In  the  equations  by  difference  operators 
and  satisfying  the  governing  equations  at  certain  discrete  points,  as  In 
the  finite  difference  method,  or  dividing  the  continuous  domain  into 
subdomains,  each  having  a  finite  number  of  degree  of  freedom,  as  In  the 
finite  element  method. 

The  boundary  element  method  (BEM), sometimes  called  the  boundary 
integral  equation  method  (BIEM),  follows  a  different  approach.  Instead 
of  working  directly  on  the  original  set  of  differential  equations,  these 
governing  differential  equations  are  converted  Into  a  set  of  Integral 
equations  before  applying  any  discretization  or  approximation  scheme. 
The  resulting  Integral  equations  to  be  discretized  Involve  only  the  val¬ 
ue  of  the  variables  at  the  extreme  range  of  the  Integration,  i.e.  on  the 
boundaries  of  the  region.  This  Implies  that  the  subsequent  discretiza¬ 
tion  scheme  needs  to  be  done  only  at  the  boundary  of  the  domain  and  the 
spatial  dimensions  of  the  problem  are  reduced  by  one.  The  distribution 
function  for  the  field  variables  over  each  discretized  element  can  be 
uniform,  linear,  or  higher  order,  similarly  to  the  finite  element  meth¬ 
od. 

The  boundary  element  method  resulted  from  the  advancements  made 
In  the  theory  of  Integral  equations  and  the  development  of  digital  com¬ 
puters.  In  the  pre-computer  era,  the  study  of  Integral  equations  and 
their  application  to  various  types  of  boundary  and  Initial  value  prob- 
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lems  had  been  a  topic  of  Interest  to  mathematicians  for  years.  The  first 
rigorous  Investigations  on  the  general  theory  of  the  subject  was  made  by 
Volterra  (1860-1940)  and  later  by  Fredholm  (1866-1927)  who  also  treated 
the  subject  In  conjunction  with  the  theory  of  linear  elastostatlcs. 
Since  then  Integral  equations  have  been  studied  Intensively  by  mathe¬ 
maticians  and  Investigators  In  the  field  of  theoretical  mechanics,  nota¬ 
bly  among  Russian  workers  as  Muskhellshvlll ,  Mlkhlln  and  Kupradze. 
Recent  analytical  work  along  these  lines  has  been  done  In  the  areas  of 
steady  state  elastodynamlcs  (Kupradze,  1963)  and  transient  elastodynam- 
Ics  (Doyle,  1966).  With  the  development  of  high  speed  digital  computers, 
many  Investigations  have  been  carried  out  to  obtain  boundary  element 
formulations  and  to  Implement  a  general  numerical  algorithm  for  solving 
a  wide  range  of  practical  problems. 

Boundary  element  methods  are  classified  by  most  Investigators 
into  three  different  types: 

1.  Indirect  Formulation.  In  this  formulation,  the  integral  equations 
are  formulated  in  terms  of  a  unit  singular  solution  that  satisfies 
the  original  differential  equation  such  as  the  free-space  (infinite 
domain)  Green's  function.  This  unit  singular  solution  Is  distrib¬ 
uted  over  the  boundary  of  the  region  In  a  manner  specified  by  density 
functions  which  are  unknowns  and  have  to  be  obtained  in  the  numerical 
solution  process.  Once  the  density  functions  are  known,  the  sol¬ 
ution  anywhere  within  the  body  can  be  calculated  by  a  simple  inte¬ 
gration  process.  Algorithms  based  on  the  Indirect  formulation  were 


presented  In  Massonet  (1965)  and  Oliveira  (1968).  Banerjee  (1971); 
and,  Tomlin  and  Butterfield  (1974)  extended  and  Improved  the  algo¬ 
rithm  to  two  and  three  dimensional  problems  Involving  burled  founda¬ 
tion  as  well  as  nonhomogeneity  and  orthotropy  of  the  medium. 

Direct  Formulation.  The  unknown  functions  In  the  Integral  equations 
of  this  formulation  are  the  actual  physical  variables  of  the  prob¬ 
lem.  The  algorithm  was  first  developed  by  Rizzo  (1967)  for  the  elas- 
tostatlc  problem;  Cruse  (1968),  Cruse  and  Rizzo  (1968)  applied  It  to 
transient  elastodynamlc  problems.  Using  Bettl-Maxwel 1  reciprocal 
theorem,  they  derived  Somlgllana's  Identity  for  the  displacement 
Inside  the  body  due  to  known  surface  displacements  and  tractions. 
The  Integral  equation  relating  displacements  and  tractions  on  the 
boundary  was  then  obtained  by  taking  the  limiting  case  of 
Somlgllana's  Identity.  This  direct  formulation  can  also  be  derived 
as  a  special  case  of  a  generalized  FEM  based  on  the  Galerkln  weighted 
residual  technique  as  In  Zlenklewlcz  (1978),  Brebbla  (1978),  or 
Brebbla  and  Walker  (1980). 

Somi-diroct  Formulation.  In  this  formulation,  the  Integral 
equations  are  formulated  In  terms  of  selected  unknown  functions  such 
as  stress  functions  or  stream  functions  for  elasticity  and  flow 
problem  respectively.  The  physical  solution  can  then  be  obtained  from 
these  unknown  functions  by  simple  mathematical  operation.  This  for- 
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initiation  was  used  In  the  work  of  Jaswon  and  Ponter  (1963)  and  Symm 
(1963). 

The  boundary  element  method  has  become  very  popular  and  has  been 
successfully  applied  to  a  variety  of  linear  problems  In  elastostatlcs, 
elastodynamlcs,  potential  theory,  fracture  mechanics,  fluid-structure 
Interaction,  soil  structure  Interaction,  porous  media  flow,  etc.  A  com¬ 
prehensive  account  of  these  applications  was  collected  In  the  books  by 
Cruse  and  Rizzo  (1975),  Brebbla  (1978),  Brebbla  and  Walker  (1980),  Breb- 
bia  (1981),  Banerjee  and  Butterfield  (1981)  and  Liggett  and  Liu  (1983). 

For  the  transient  elastodynamlc  problem,  the  method  Is  very  use¬ 
ful  especially  In  the  field  of  geophysics  and  dynamic  soil  structure 
Interaction.  In  this  case,  the  boundary  element  method  can  be  applied  In 
conjunction  either  with  Integration  In  the  time  domain  as  In  Cole  et  al 
(1978),  or  with  frequency  domain  analysis  using  Laplace  or  Fourier 
transform  as  In  Dominguez  (1978).  With  the  recent  development  of  the 
dynamic  Green's  function  for  a  half  space  (Johnson,  1974),  and  a  layered 
half  space  (Luco  and  Apsel,  1983),  the  research  on  the  application  of  the 
boundary  element  method  In  dynamic  geotechnical  engineering  and  geophy¬ 
sics  appears  to  be  a  very  promising  area. 
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CHAPTER  3 


FINITE  ELEMENT  FORMULATION 


To  solve  the  problem  of  wave  propagation  In  an  infinite  medium 
with  Inclusions  using  finite  elements,  It  Is  necessary  to  represent  an 
Infinite  domain  with  a  finite  one.  This  approximation  can  give  good 
results  If  the  boundaries  of  the  domain  are  extended  far  enough  so  that 
the  amplitude  of  the  reflected  waves  reaching  the  region  of  Interest  Is 
very  small  (If  there  Is  Internal  damping),  or  so  that  the  time  of  arrival 
of  these  reflected  waves  exceeds  the  duration  of  Interest.  At  the  same 
time,  the  size  of  the  finite  element  mesh  should  be  fine  enough  to 
reproduce  accurately  the  propagation  of  waves  of  a  given  wavelength. 
With  these  constraints.  If  a  formulation  that  requires  the  assembly  of 
the  global  stiffness  matrix  Is  used.  It  will  result  In  an  extremely  large 
system  of  equations  which  Is  Impractical  to  handle.  With  some  simplified 
assumptions  Involving  the  use  of  a  diagonal  mass  matrix  and  neglecting 
velocity-dependent  damping  forces,  It  Is  possible,  however,  to  solve  the 
problem  using  an  explicit  Integration  scheme,  such  as  the  central  dif¬ 
ference  method,  and  then  marching  out  the  solution  In  the  time  domain. 
The  solution  can  essentially  be  carried  out  at  the  element  level  without 
assembling  the  global  stiffness  matrix.  This  numerical  scheme  can  dras¬ 
tically  reduce  the  size  of  computer  storage  required.  The  solution  Is, 
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still,  time  consuming  since  the  time  step  of  Integration  must  be  kept 
sufficiently  small  to  guarantee  stability. 

To  keep  the  computation  cost  within  a  reasonable  limit,  only  a 
two-dimensional  finite  element  model  with  four-node  linear  elements  will 
be  used  and  presented  In  this  chapter. 


3.1  FINITE  ELEMENT  MODEL 


The  concept  of  the  finite  element  method  Is  very  broad  and  can  be 
applied  In  a  number  of  different  formulations.  In  this  work,  the  stiff¬ 
ness  or  displacement  method  Is  used.  This  method  Is  based  on  assuming 
functions  to  approximate  the  displacement  field  for  each  element.  These 
displacement  functions  can  be  expressed  In  various  forms,  the  simplest 
linear  functions  for  four-node  rectangular  element  will  be  used  here. 
The  equations  of  equilibrium  for  a  finite  element  system  In  motion  can 
then  be  derived  either  by  using  the  principle  of  virtual  displacement  In 
conjunction  with  d'Alembert  forces  or  by  applying  Hamilton's  principle. 


3.2  EQUATION  OF  EQUILIBRIUM 

With  the  viscous  damping  term  neglected  for  simplicity,  the 
equation  of  equilibrium  derived  for  a  finite  element  system  can  be  writ¬ 


ten  as 
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MU  +  KU  =  R 


(3.1) 


where  M  and  K  are  mass  and  stiffness  matrices;  R  Is  the  external  load 
vector;  U  and  U  are  the  displacement  and  acceleration  vectors  of  the 
finite  element  assemblage. 

The  stiffness  matrix  K  Is  formed  by  assembling  each  element 
e 

stiffness  matrix  K  which  Is  expressed  as 

Ke  ■  f  BTDB  dfte  (3.2) 

where  B  Is  a  matrix  relating  element  strains  to  nodal  displacements;  D  Is 
a  matrix  relating  stresses  to  strains  In  the  element  and  sometimes  called 
the  element  elasticity  matrix.  The  Integral  above  Is  carried  out  over 

p 

the  region  Q  of  the  element. 

The  only  external  load  considered  In  our  case  Is  the  concen¬ 
trated  load  at  the  excitation  point.  The  element  mass  matrix  Me  Is  In 
the  form 


Me  *  [  PNTN  dne 


(3.3) 


where  p  Is  the  mass  density,  and  N  denotes  the  matrix  of  assumed  dls 
placement  functions. 

The  element  mass  matrix  resulting  from  (3.3)  is  called  the  con 


si  stent  mass  matrix  because  the  same  displacement  functions  (sometimes, 
called  the  basis  functions)  are  used  to  formulate  the  stiffness  and  mass 
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matrices.  If  this  consistent  mass  matrix  Is  used,  the  resulting  global 
mass  matrix  M  will  be  non-diagonal.  It  Is,  however,  not  a  requirement 
for  convergence  to  use  a  consistent  mass  matrix;  therefore,  different 
functions  may  also  be  used.  If  the  basis  functions  that  have  unit  value 
over  the  region  tributary  to  the  node  under  consideration  and  zero  else¬ 
where  are  used  and  the  tributary  regions  do  not  overlap,  the  resulting 
matrix  M  will  be  diagonal.  The  use  of  a  diagonal  (or  lumped)  mass  matrix 
Involves  much  less  computational  effort  and  suffers  no  loss  In  order  of 
convergence  and  accuracy  In  most  cases  especially  In  wave  propagation 
problems  with  linear  displacement  expansions  (Scalettl,  1977). 

For  the  finite  element  numerical  scheme  used  In  this  work.,  It  Is 
necessary  to  use  a  diagonal  mass  matrix  because  the  solutions  will  be 
carried  out  on  the  element  level  without  assembling  the  global  matrices. 
Furthermore,  the  benefits  of  a  diagonal  mass  matrix  are  evident  when 
applying  a  conditionally  stable  Integration  scheme  since  the  element 
will  have  a  lower  natural  frequency  (l.e.,  larger  natural  period)  and  a 
larger  time  step  of  Integration  may  be  used. 


3.3  DIRECT  INTEGRATION  METHOD 


Mathematically,  Equation  (3.1)  represents  a  system  of  linear 
differential  equations  of  second  order  which  can  be  solved  by  Inte¬ 
gration  using  a  numerical  step-by-step  procedure.  The  numerical  Inte¬ 
gration  Is  based  on  approximations.  First,  the  equilibrium  equation 
(3.1)  Is  satisfied  only  at  each  discrete  time  (with  Interval  At)  Instead 


of  at  any  time  t;  second,  a  variation  of  displacements,  velocities,  and 
accelerations  within  each  time  Interval  At  Is  assumed.  The  first  one  Is 
the  basis  for  all  finite  difference  methods.  The  latter  concept  differ¬ 
entiates  one  method  from  another  (such  as  central  difference  and  forward 
difference  methods,  etc.)  and  determines  the  accuracy,  stability,  and 
cost  of  the  solution  procedure. 

One  of  the  procedures  that  Is  effective  and  will  be  used  In  this 
work  Is  the  central  difference  method  which  assumes  that 


<  Vt-  2V 


IF" 


ut-u> 


(3.4) 


where  the  subscripts  In  the  above  equation  denote  the  discrete  time  at 
which  the  values  are  obtained. 

The  velocity  expansion  which  has  the  same  order  of  error  Is 


*Ut+At~  Ut-At* 
t  *  2At 


(3.5) 


Considering  the  equilibrium  equation  (3.1)  at  time  t 


MUt  +  KUt  -  Rt  (3.6) 

and  substituting  the  relations  for  In  (3.4)  Into  (3.6),  we  get 

<ZFM>lW  *  Rt  "  (K  " 


(3.7) 
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The  method  assumes  that  the  Initial  conditions  (the  displace- 
ment,  velocity,  and  acceleration  at  time  t  =  0)  are  known.  The  solutions 
at  the  next  time  steps  are  then  calculated  using  Equation  (3.7)  up  to  the 
end  of  the  time  span  under  consideration.  Since  the  solution  of  U. ...  Is 

X+aX 

calculated  using  the  equilibrium  conditions  at  time  t  Instead  of  at  time 
t+At,  the  method  Is  called  an  explicit  Integration  scheme. 

Equation  (3.7)  can  be  rewritten  In  the  form  of  summations  of  the 
element  matrices 

l  ak"6  Ut*At  ■  -  1(6  Ut>3  ♦  Cl  (2Ut  -  Ut.at)]  (3.8) 

6  6  6 

If  the  global  mass  matrix  M  is  diagonal,  Equation  (3.7)  can  be 
solved  without  factoring  a  matrix,  and  only  matrix  multiplication  are 
necessary  on  the  right  hand  side  of  (3.7).  Let  m^  denotes  the  ith  diag¬ 
onal  element  of  the  mass  matrix  M,  the  1th  components  of  the  vectors 
Ut+At  can  be  determined  from  (3.8)  as 

uLt  ■  -  kV’  *  2ui  -  uLt  cj.s) 

6  6 

where  the  summation  term  Indicates  the  result  of  R^-K  Ut  for  the  1th  com¬ 
ponents  contributed  by  each  element  'e'.  Therefore,  all  calculations 
can  be  done  at  the  element  level  and  the  solutions  can  be  marched  out  In 
the  time  domain  without  assembling  the  global  stiffness  matrix.  The 
method  becomes  even  more  effective  when  the  element  stiffness  matrices 
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K  of  all  elements  are  the  same  since  only  one  element  stiffness  matrix 
needs  to  be  calculated.  It  can  be  stored  and  then  recalled  for  uses. 

In  the  formulation  shown  In  (3.9),  It  Is  possible  to  include  the 
damping  terms  without  a  significant  Increase  In  the  cost  of  computation 
If  the  damping  matrix  Is  diagonal.  With  this  condition,  the  benefits  of 
performing  the  solution  on  the  element  level  are  still  preserved.  Since 
the  main  purpose  of  the  finite  element  formulation  is  to  check  and  com¬ 
pare  with  the  boundary  element  method  (BEM)  and  since  the  same  type  of 
damping  cannot  be  Incorporated  Into  both  models,  the  damping  terms  are, 
however,  neglected. 


3.4  TIME  STEP  OF  INTEGRATION 

The  central  difference  method  used  In  the  formulation  presented 
In  Section  3.3  Is  classified  as  a  conditionally  stable  integration 
scheme.  It  requires  that  the  time  step  At  be  smaller  than  a  critical 
value  Atcr  to  guarantee  stability,  l.e 

At  s  At  =  l  (3.10) 
cr  it 

where  T  Is  the  minimum  period  of  the  finite  element  assemblage.  This 
minimum  period  T  can  be  calculated  using  any  of  the  numerical  techniques 
available  for  the  solution  of  elgenproblems,  which  Is  not  possible  In  our 
case.  Since  the  solution  for  a  large  eigenvalue  problem  Is  also  time 
consuming  and  expensive,  It  is  not  practical  to  solve  the  large  elgenva- 
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lue  problem  just  to  determine  the  appropriate  time  step  At  for  the  direct 
integration  method.  A  simpler  method  Is  to  determine  a  lower  bound  for  T 
which  Is  controlled  by  the  smallest  natural  period  of  any  element  in  the 
finite  element  model.  The  minimum  period  for  a  single  element  can  be 
calculated  numerically  from  the  stiffness  and  mass  matrices  of  the  ele¬ 
ment  or  It  can  be  determined  from 

Tmin  ‘  c,c2/(6.5  -  v)  ±  (3.11) 

The  equation  above  applies  for  a  four-node,  rectangular,  plane 
strain  element  where  t  Is  the  minimum  element  dimension,  c$  is  the  shear 
wave  velocity,  v  Is  Poisson's  ratio,  Cj  is  either  it  for  lumped  masses  or 
v/S3  for  consistent  masses,  and  is  a  shape  factor,  equals  to  1  for  a 
square  element. 

The  benefit  of  using  lumped  masses  over  a  consistent  mass  matrix 
mentioned  above  Is  evidenced  by  Equation  (3.11).  The  lumped  mass  approx¬ 
imation  results  In  a  larger  natural  period,  hence  a  larger  time  step  At 
may  be  used.  On  the  other  hand,  any  reduction  In  the  size  of  one  single 
element  (which  will  reduce  a  certain  diagonal  element  In  a  global  mass 
matrix)  will  decrease  the  minimum  period  T  of  the  system.  Any  increases 
In  the  stiffness  (or  shear  wave  velocity  c$)  on  one  element  will  have  the 
same  effect. 


3.5  FINITE  ELEMENT  MESH  AND  BOUNDARY  CONDITIONS 


In  this  work,  four-node  square  elements  with  linear  basis  func¬ 
tions  will  be  used.  All  elements  used  In  the  model  are  the  same  in  order 
to  take  advantage  of  the  solution  scheme  mentioned  earlier  when  all  ele- 
ments  have  the  same  stiffness  matrices  K  .  Each  Inclusion  is  represented 
by  a  single  element  and  Is  assumed  to  be  perfectly  rigid  with  compatibil¬ 
ity  satisfied  everywhere.  It  Is  also  possible  to  modify  the  model  for 
inclusions  of  varying  stiffness  but  the  time  step  of  integration  must  be 
changed  according  to  the  stability  crlterlor  mentioned  In  Section  3.4. 

The  size  of  a  finite  element  mesh  relative  to  the  wavelength 
determines  the  accuracy  of  the  results.  Typically,  nodal  spacing  In  the 
direction  of  the  propagating  wave  should  range  from  1/4  to  1/10  of  the 
shortest  wavelength  of  Interest,  and  the  value  of  1/8  will  be  used 
throughout  this  work. 

For  a  finite  element  model,  an  Infinite  medium  Is  approximated 
by  a  finite  domain  with  fixed  or  free  boundaries.  These  boundaries  may 
hamper  any  motions  In  their  vicinity;  In  addition,  the  reflected  waves 
from  these  boundaries  will  be  combined  with  the  direct  wave  making  inter¬ 
pretation  of  the  results  more  difficult.  Therefore,  the  domain  has  to  be 
sufficiently  large  to  make  boundary  effects  negligible  over  the  time 
span  of  Interest.  The  Investigation  for  the  proper  size  of  the  domain 
will  be  discussed  In  Chapter  5. 

Since  each  Inclusion  Is  to  be  represented  by  a  single  element, 
the  centerline  of  the  domain,  which  Is  also  the  principal  axis  of  wave 


propagation,  has  to  pass  through  the  center  of  Inclusions  (or  element)  In 
order  to  take  advantages  of  symmetry  and  antisymmetry.  Because  of  these 
conditions,  the  excitation  at  a  point  has  to  be  approximated  by  distrib¬ 
uting  the  force  evenly  between  two  nodes  as  shown  In  the  model  of  Figure 
3.1.  This  approximation  Is  not  serious  as  long  as  the  distance  between 
these  two  nodes  Is  small  compared  with  the  distances  from  the  excitation 
point  to  the  Inclusion  and  the  target  point. 


CHAPTER  4 


BOUNDARY  ELEMENT  FORMULATION 

The  boundary  element  method  Is  basically  a  solution  In  discrete 
form  of  the  boundary  Integral  equation  that  relates  the  tractions  and  the 
displacements  over  the  boundaries  of  the  domain  under  consideration. 
The  discretization  Is  achieved  by  dividing  the  boundaries  Into  a  number 
of  smaller  elements  and  then  calculating  the  Integrals  numerically.  The 
variations  of  the  displacements  and  the  tractions  may  be  assumed  to  be 
constant  or  of  a  higher  order  (  linear,  quadratic,  cubic,  etc.).  This 
procedure  leads  to  a  system  of  equations  with  the  tractions,  the  dis¬ 
placements  or  both  as  unknowns.  Once  these  tractions  and  displacements 
over  the  boundary  are  solved  for,  the  results  at  any  point  Inside  the 
domain  can  be  calculated. 

For  time-dependent  problems,  the  boundary  element  method  may  be 
solved  either  In  the  time  domain  by  a  step-by-step  Integration  (Cole  et 
al,  1978)  or  In  the  frequency  domain  by  eliminating  the  time  variable 
with  a  Laplace  transform  (Cruse  and  Rizzo,  1968)  or  a  Fourier  transform 
(Dominguez,  1978).  Once  the  equations  are  solved  In  the  transformed 
space,  the  original  variable  Involving  the  time  dependence  may  be  recov¬ 
ered  by  Inverting  the  transform  numerically. 
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The  boundary  element  formulation  described  In  this  chapter  Is 
based  on  the  method  used  by  Dominguez  (1978),  who  studied  the  dynamic 
stiffness  of  surface  and  embedded  foundations.  For  the  case  of  wave  pro¬ 
pagation  through  an  Infinite  medium  with  Inclusions,  some  modifications 
are  needed  In  the  formulation.  By  working  In  the  frequency  domain  using 
the  Betti-Maxwell  reciprocity  relation  and  a  fundamental  solution  In 
transform  space  for  a  unit  harmonic  point  force  In  an  Infinite  medium,  an 
equation  for  displacement  Inside  the  body  due  to  the  known  boundary  trac¬ 
tions  and  displacements  (so  called  Somlgllana's  Identity)  Is  derived. 
The  integral  equation  Is  then  formed  by  taking  the  limiting  case  of 
Somlgllana's  Identity  for  points  located  of  the  boundary  of  the  domain. 
By  discretizing  the  boundaries  Into  small  elements,  the  boundary  matri¬ 
ces  representing  a  system  of  equations  that  relate  the  displacements  and 
the  tractions  at  the  boundaries  are  formed.  This  formulation  Is  carried 
out  In  two  separate  parts:  (1)  for  the  Inclusions;  and  (2)  for  the  infi¬ 
nite  medium  with  cavities  and  an  excitation  by  a  unit  harmonic  point 
force.  These  two  parts  are  then  combined  using  equilibrium  and  compat¬ 
ibility  conditions.  This  leads  to  a  system  of  equations  which  can  be 
solved  for  unknown  displacements  and  tractions  at  the  surface  of  the 
Inclusions.  Then,  the  displacements  at  any  points  Inside  the  medium  are 
calculated  using  Somlgllana's  Identity.  For  a  frequency  domain  analy¬ 
sis,  these  displacements  or  transfer  functions  are  solved  for  different 
frequencies.  The  solution  for  a  transient  excitation  can  then  be 
obtained  using  Fourier  transforms. 


.■  y  1 
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The  details  of  the  boundary  element  formulation  are  presented  In 
the  following  sections. 


4.1  BOUNDARY  ELEMENT  MODELS 


4.1.1  Model  for  Soil  with  Cavity 


The  formulation  starts  with  the  Betti-Maxwell  reciprocity 
relation  between  the  Fourier  transforms  of  the  actual  and  virtual 
states.  Consider  a  body  Q  with  a  boundary  r  under  a  certain  system  of 
loads  (see  Figure  4.1a).  Assuming  a  virtual  state  on  the  same  body, 
the  following  relation  can  be  established 


f  -  r 

A  <N* 

*4* 

JrVi dr  +  J, 

b.u.  da 

i  ’  1 

LVi dr  *  J 

aVi 

da 


(4.1) 


The  usual  Indlclal  notation  of  Cartesian  tensor  analysis  Is  used. 

The  variables  u^,  t^,  and  b^  are  the  components  of  displacement,  traction 

and  body  force,  respectively.  The  sign  1  ■*'  over  each  terms  Indicates 

that  they  are  In  the  transform  space  (l.e.,  the  frequency  domain). 

Let  each  component  of  the  virtual  force  b^  be  a  Dirac  delta  f unc¬ 
le 

tlon  A  which  represents  a  unit  concentrated  force  at  point  k,  then 


a* 

b1 


A  e. 


(4.2) 
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(a)  Infinite  Medium  with  an  Interior  Boundary  (Cavity)  and 
an  Excitation  Force  as  a  Body  Force 


Boundary  (Inclusion) 


Figure  4.1  -  Virtual  State  for  an  Infinite  Medium 
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where  e^  Is  a  unit  normal  vector  In  the  1  direction. 

The  virtual  state,  therefore,  becomes  the  fundamental  solution, 
or  the  response  of  the  Infinite  medium  to  a  unit  concentrated  harmonic 
force  applied  at  point  P.  This  solution,  following  the  works  by  Doyle 
(1966)  and  Cruse  and  Rizzo  (1968),  can  be  obtained  from  the  equations  of 
motion  of  a  linear  Isotropic  elastic  body  and  Is  of  the  form 


“ji  -skf  K-j  -  xrMr>j] 

Tji  *  sr  c<&  -  &  ♦  '■•i-j)  -  -  2rMr-j  a? 


(4.3) 


(4.4) 


A  A 

where  and  T^  are  respectively  the  displacement  and  traction  In  1 
direction  at  point  i  due  to  the  unit  force  In  direction  j  at  point  k,  r  Is 
the  distance  from  k  to  t,  p  Is  the  mass  density,  c$  and  cp  are  the  veloci¬ 
ty  of  shear  wave  and  compression  wave  of  the  medium  respectively,  and  n 
Is  a  unit  outer  normal  vector. 

In  the  equations  above,  a  =  2  for  the  two-dimensional  case,  and  a 
=  4  for  three  dimensions.  The  functions  ♦  and  x  for  the 
three-dimensional  case  are 


*  * 


x  = 


e(-1o)r/cs) 

r 

e(-1a»r/c$) 

_ 


(-1u>r/cp) 

r 

e(-1u>r/cp) 

r 


(4.5) 

(4.6) 
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For  the  two-dimensional  case,  i>  and  x  are 


*■  +  -d 


—  u  /  j \  S  1/  /  "lliif*  \ 
X  =  k2(— )  -  ^  k2(~) 

s  p  p 


P  P 


(4.7) 

(4.8) 


where  Kg,  K^,  and  <2  are  modified  Bessel  functions  of  the  second  kind  of 
orders  zero,  one,  and  two  respectively.  Note  that 


ui  =  Uj1  ej 


*1  *  Tji  ej 


(4.9) 


by  substituting  (4.9)  and  (4.2)  Into  (4.1),  we  get 


fr*iVj dr  *  d!1  ’  dr  *  IQVk  ei d0  (410) 


The  dummy  index  of  the  last  term  of  (4.10)  can  be  changed  from  1 


to  j,  and 


'n 


A  k 

UjA  dn 


(4.11) 


k  A  k 

where  A  and  uj  are  the  delta  function  and  displacement  at  point  k  In 
direction  j  respectively.  Equation  (4.10)  becomes 
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"5  ■  dr  -  frvi dr  *  6a  ‘4-i2> 

Equation  (4.12)  Is  known  as  Somlgllana's  Identity,  which  will  be 
used  later  to  compute  the  displacement  at  any  Interior  point.  In  order 
to  formulate  the  problem  as  a  boundary  technique,  the  virtual  load  point 
k  has  to  be  taken  to  the  boundary.  Consequently,  the  Integral  along  the 
boundary  will  present  a  singular  point.  This  singular  Integration  can  be 
done  by  considering  the  hemisphere  on  the  boundary  of  a 
three-dimensional  model  (semi-circle  for  the  two-dimensional  case).  The 
boundary  point  k  Is  assumed  to  be  at  the  center  of  the  sphere,  and  the 
radius  'e'  Is  reduced  to  zero.  The  point  will  then  become  a  boundary 
point. 

For  a  smooth  boundary  at  the  point  k,  It  can  be  shown  that 


and 


11m 

e  -*■  0 


11m 

e  -*■  0 


lrVt dr 

e 

lrVi dr 

e 


=  0 


where  e  Is  the  radius  of  the  hemisphere  over  point  k.  The  Somlgllana's 
Identity  (4.12)  can  then  be  written  for  a  boundary  point  as 


.Vi 


dr 


dn 


(4.13) 
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Equation  (4.13)  will  now  be  applied  on  the  boundary  of  the  domain 
under  consideration.  For  the  wave  propagation  model,  there  are  no  body 
forces  except  at  the  excitation  point,  where  b^  Is  a  unit  Impulse  func¬ 
tion.  Therefore,  the  last  term  of  (4.13)  becomes 


j^JI  dD  -  <“o>j  «-14> 

a  k 

where  (uQ)j  Is  theoretically  the  displacement  at  the  excitation  point 
(In  the  same  direction  as  the  excitation  force)  due  to  the  virtual  load 
at  point  k  In  j  direction.  Since  both  the  excitation  force  and  the  vir¬ 
tual  load  are  the  same  (unit  point  load),  by  the  reciprocal  theorem, 

A 

(uQ)j.  Is  also  equal  to  the  displacement  at  the  point  k  in  the  j  direction 
due  to  the  excitation  force. 

The  solution  of  the  Integral  equation  (4.13)  can  be  done  In  a 
systematic  way  using  numerical  procedures  by  dividing  the  boundary  T 
Into  n  elements.  In  this  work,  the  constant  element  Is  used  (l.e.,  there 
Is  only  one  node  at  the  center  of  each  element).  The  displacements  and 
tractions  are  assumed  to  be  constant  throughout  each  element  and  equal  to 
the  value  at  the  mid-node. 

By  considering  the  discretized  model  as  shown  In  Figure  4.2,  and 
using  Equation  (4.14),  Equation  (4.13)  can  be  written  In  discretized 
form  for  a  given  point  k  as 


dr  {y, 


dr  <y, 


<%>$ 


(4.15) 
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k:  Nodal  point  under  consideration  where  virtual  load  is  applied 
a:  Element  over  which  the  integration  is  carried  out. 


Figure  4.2  -  Discretized  Model 


t 


The  above  equation  applies  for  a  particular  node  k.  Define  the 


terms 


Ujl  dr 


^kz^ji 


(4.16) 


Both  Integrals  above 
the  Integral  Is  carried  out. 


=  ^kz]ji 

relate  the  k  node  with  element  Z  over  which 
Equation  (4.15)  can  be  rewritten  as 


1~k  .  , 

*  l 

d  J  Z=1 


n 

l  [H 


kz-Jji{uz}i 


r> 

I  CG, 

Z=1 


kz^j7{tz^i 


(4.17) 


Equation  (4.17)  relates  the  value  of  u  at  mid-node  k  with  the 
value  of  u  and  t  at  all  the  nodes  on  the  boundary,  Including  k.  It  can  be 
written  for  each  k  node  obtaining  n  equations.  Let 


Hkl  =■  Hm  -hen 

Hia  ’  Sk«*  1  when  k  '  1 

Equation  (4.17)  can  then  be  written  as 


(4.18) 


£-CHkz]j1{Vi 


z=l 


i  &  + 


(4.19) 


The  complete  system  of  equations  for  the  n  nodes  can  be  expressed 


In  matrix  form  as 
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HU  =  GT  +  UQ  (4.20) 

where  U  and  T  are  vectors  with  two  (three  for  3-D  case)  displacement  and 
traction  components  for  each  of  the  n  nodes.  UQ  is  the  vector  of  free 
field  displacement  components  at  each  of  the  n  nodes  due  to  the  unit 
excitation  force  and  can  be  calculated  directly  from  Equation  (4.3).  G 
and  H  are  2nx2n  square  matrices  (3n*3n  for  3-D  case)  obtained  by  assembl¬ 
ing  the  2*2  submatrices  (3x3  for  3-D  case)  that  relate  every  two  element 
nodes  k  and  l 

Note  that  the  matrices  G  and  H  are  not  symmetric  and  most  of 
their  elements  differ  from  zero.  In  spite  of  these  factors,  the  boundary 
element  method  is,  in  many  cases,  very  efficient  due  to  the  much  smaller 
number  of  unknowns.  This  is  because  the  spatial  dimension  of  the  problem 
is  reduced,  and  only  the  boundary  of  the  region  under  study  needs  to  be 
di scretlzed. 


4.1.2  Model  for  Inclusion 

Equation  (4.20)  Is  formulated  as  a  relationship  between  the  dis¬ 
placements  and  the  tractions  at  the  boundary  of  the  cavity  In  an  Infinite 
medium  with  a  unit  excitation  force.  In  order  to  formulate  the  problems 
of  wave  propagation  In  an  Infinite  medium  with  Inclusions,  a  similar 
relationship  Is  needed  for  each  Inclusion.  The  equations  for  the  Infi¬ 
nite  medium  and  the  Inclusions  are  then  combined  from  the  equilibrium 
condition  and  the  compatibility  of  displacements  at  Interfaces. 


36 


The  derivation  of  an  equation  like  (4.20)  for  an  Inclusion  fol¬ 
lows  the  same  procedure  as  In  Section  4.1.1  except  that  there  are  no  body 
force  terms  In  this  case.  By  considering  the  model  shown  In  Figure  4.1b, 
starting  with  Equation  (4.1),  and  neglecting  the  body  force  term  In  the 
actual  state,  a  system  of  equations  that  establishes  a  relation  between 
displacements  and  tractions  over  the  surface  of  an  Inclusion  Is  obtained 
as 


Vb  ■  Vb  <4-21> 

All  symbols  In  the  equation  above  are  as  defined  for  Equation 
(4.20),  while  the  subscript  ' b'  Indicates  that  they  belong  to  an  Inclu¬ 
sion. 

It  Is  worth  noting  here  the  distinctive  structure  of  the  bounda¬ 
ry  matrices  (G  and  H  matrices)  for  the  cavity  and  for  the  Inclusion.  In 
the  boundary  element  formulation,  the  major  difference  between  an  inclu¬ 
sion  and  a  cavity  Is  the  direction  of  the  outer  normal  vector  n  used  In 
the  calculation  of  traction  as  shown  In  Equation  (4.4).  It  can  be 
seen  from  Equation  (4.16)  that  the  calculation  of  the  G  matrix  Involves 

A 

the  Integration  of  the  displacement  which  Is  Independent  of  n;  there¬ 
fore,  the  matrix  G  (for  soil  with  a  cavity)  and  the  matrix  G^  (for  an 
Inclusion)  are  Identical  provided  that  their  surfaces,  their  elements, 
and  their  properties  are  the  same.  In  the  same  way,  the  H  matrix  Is 
obtained  by  Integrating  the  traction  tji  which  depends  on  the  direction 
of  the  outer  normal  n.  Consequently,  the  matrix  H  (for  soil  with  a  cavl- 


ty)  and  the  matrix  Hb  (for  an  Inclusion)  will  have  the  same  value  but 
opposite  sign  except  for  each  diagonal  elements  of  the  matrix  which  are 
still  the  same  and,  for  our  case,  equals  to  1/2.  In  summary, 

G  =  Gb 

and 

H  =  I  -  Hb 

where  I  Is  an  Identity  matrix. 


4.1.3  Modal  for  Wava  Propagation  in  a  Madium  with  inclusions 


A  formulation  for  wave  propagation  In  an  Infinite  medium  elth 
Inclusions  Is  derived  by  combining  (4.20)  and  (4.21)  under  two  condi¬ 
tions  Imposed  at  the  Interfaces  of  soil  and  Inclusions.  These  two  condi¬ 
tions  are,  first,  that  the  tractions  at  Interfaces  must  be  In 
equilibrium,  and  second,  that  there  Is  no  separation  between  soil  and 
Inclusions  at  those  Interfaces.  Both  conditions  Imply  that 


and 


Tb  =  -T 
ub  =  U 


(4.22) 


Equation  (4.21)  can  be  rewritten  as 
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By  using  (4.22),  Equation  (4.23)  becomes 

T  =  -G^H^  (4.24) 

Substituting  (4.24)  Into  (4.20)  and  rearranging,  an  equation  Is 
obtained  that  can  be  solved  for  the  displacements  at  the  surfaces  of  the 
Inclusions 


(H*GG-\)U  =  U0 


(4.25) 


In  the  above  equation,  Gj^H^  belongs  to  each  Inclusion  and  has 
the  same  size  as  G  and  H  If  there  Is  only  one  Inclusion  In  the  system. 
For  the  case  of  more  than  one  Inclusion,  the  matrix  Gj^Hb  for  each  Inclu¬ 
sion  will  be  a  submatrix  along  the  diagonal  of  the  global  matrix  because 
the  Inclusions  cannot  Interact  with  each  other  except  through  the  soil. 
The  global  matrix  Gj^Hb  for  the  case  of  multiple  Inclusions  Is  typified 
In  Figure  4.3. 

Once  the  displacements  U  at  boundary  are  solved,  the  boundary 
tractions  T  can  be  calculated  using  Equation  (4.24).  The  displacement  at 
any  point  In  the  domain  can  then  be  computed  using  equation  (4.12),  which 
can  be  written  for  the  displacement  at  any  Interior  point  k  In  dlscrea- 
tlzed  form  as 


n 

J, 


tGk*^j1{Vl 


J,  tSwVVi  *  ("o'j 


(4.26) 
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In  the  above  equation,  and  [H^]^  are  the  Integrals  as 

defined  In  (4.16),  but  In  this  case  they  relate  the  Interior  point  k  with 

*  k 

each  mid-node  on  the  boundary  elements.  (uQ)j  Is  the  free  field  dis¬ 
placement  component  at  the  Interior  point  k  due  to  the  excitation  force, 

A  A 

which  can  be  readily  calculated  from  (4.3).  and  t^  are  known  boundary 
displacements  and  tractions  that  have  already  been  computed  from  (4.25) 
and  (4.24). 

*  k 

It  should  be  noted  that  the  final  solution  Uj  resulting  from 
(4.26)  Is  only  the  frequency  response  for  a  unit  concentrated  harmonic 
excitation.  It  Is  also  called  a  transfer  function  and  has  to  be  deter¬ 
mined  for  several  frequencies  of  excitation  In  order  to  get  the  time 
response  from  frequency  domain  analysis. 


4.2  FREQUENCY  DOMAIN  ANALYSIS 

The  frequency  domain  approach  Involves  expressing  the  applied 
loading  In  terms  of  harmonic  components,  evaluating  the  response  of  the 
system  to  each  component,  and  then  superimposing  the  harmonic  responses 
to  obtain  the  total  response. 

Any  periodic  function  p(t)  can  be  expressed  In  a  Fourier  series 
as 


where 


P(t)  -  l 

na-» 


e(1noi1t) 


(4.27) 
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c_  -  \  fTp  P(t)  e(_in“lt)  dt  (4.28) 

n  'p  J0 

In  which  Tp  Is  a  period  of  the  loading  function,  wj  Is  the  fundamental 
frequency  and  n  Is  an  Integer. 

Defining  U(w)  as  a  frequency  response  for  a  unit  concentrated 
harmonic  excitation,  or  a  transfer  function,  It  Is  possible  to  express 
the  response  of  the  system  to  each  harmonic  component  of  the  loading 
function  which  represents  each  term  of  the  Fourier  series.  Then,  from 
the  principle  of  superposition,  the  total  response  of  the  system  to  the 
forcing  function  p(t)  can  be  written  as 

u(t)  =  l  U( noj-j )  cn  e(in“lt) 
n=-« 

A  nonperiodic  loading  can  also  be  represented  by  the  Fourier 
series  by  considering  that  the  loading  function  has  a  period  Tp  that 
extends  to  Infinity.  It  is  necessary  to  reformulate  the  Fourier  series 
expression  so  that  It  extends  over  an  Infinite  time  range.  Let 

1  “l  _  to 

Tp  ^  =  Zw 

Hut  *  nAui  =  u>_ 

I  n 


With  the  above  notation,  the  Fourier  series  (4.27)  and  (4.28)  becomes 
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P(t)  -  If  l  C(„n) 

n=-« 

t=T_/2 

c(oin)  =  |  p(t)  dt 


(4.29) 

(4.30) 


When  Tp  -*  •  ,  the  frequency  Increment  becomes  an  Infinitesimal 
(Aw  ■*  dw),  and  the  discrete  frequencies  wn  become  a  continuous  variable 
w.  Therefore,  the  Fourier  series  (4.29)  becomes  the  following  Fourier 
integral 


p(t)  =  |  c(u>)  dw 

In  which  the  harmonic  amplitude  function  Is  given  by 


(4.31) 


c(u>)  =  |°°  p(t)  eHwt)  dt  (4.32) 
t=-« 

The  two  Integrals  of  (4,31)  and  (4.32)  are  known  as  a  Fourier 
transform  pair.  Equation  (4.31)  is  the  representation  of  an  arbltary 
loading  as  an  Infinite  sum  of  harmonic  components,  where  (l/2e).c(w) 
defines  the  amplitude  per  unit  w  of  the  load  component  at  frequency  w. 
Multiplying  (4.31)  by  the  transfer  function  U(w)  yields  the  amplitude 
per  unit  w  of  the  response  components  at  frequency  w.  The  total  response 
can  then  be  obtained  by  summing  these  response  components  over  the  entire 
frequency  range,  leading  to  the  equation  for  the  analysis  of  response 
through  the  frequency  domain. 
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u(t)  =  ^  I*  H(.)  c(.)  e(wt)  d»  (4.33) 

In  order  to  make  a  practical  use  of  the  frequency  domain  method 
described  above,  It  Is  necessary  to  formulate  the  Fourier  transform  pair 
(4.31)  and  (4.32)  so  that  they  can  be  calculated  numerically.  A  discrete 
Fourier  transform  (DFT)  expression  Is,  therefore,  derived  with  the 
assumption  that  the  Input  function  Is  periodic  of  period  Tp  In  order  to 
replace  the  Infinite  time  Integral  with  the  finite  sum.  The  period  Tp 

also  serves  to  define  the  frequency  Increment  and  the  lowest  frequency 

that  may  be  considered  In  the  analysis,  thus 

w,  -  Aw  - 

P 

and  nuij  =  nAu  =  un  »  n  3  0,  1,  2 .  (4.34) 

The  loading  function  of  period  Tp  Is  then  divided  Into  N  equal 

time  increments  At,  and  the  function  Is  defined  for  each  of  the  discrete 
times  t  ,  l.e. 

Ill 

tm  =  mAt,  m  =  0,  1,  2, .  (4.35) 

With  (4.34)  and  (4.35),  the  exponential  term  In  (4.29)  becomes 

e(1«ntm)  «  e(1nAwmAt)  .  e(2ir1  p 

and  then,  (4.29)  takes  the  discrete  form 
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a  N-l  /a  j  nm \ 

p(tm)  =  §?  I  c(o.n)  e(ZlT'  N  ,  m  *  0,  1,  2 . N-l  (4.36) 

n=0 

Accordingly,  the  discrete  form  for  an  amplitude  function  c(wn)  Is 
N-l  /  n»H% 

c(o>  )  =  At  l  p(tj  e'  2  1  N  »  n  =  0,  1.  2 . N-l  (4.37) 

n  m=0  m 

Equations  (4.36)  and  (4.37)  are  the  DFT  pairs  correspond  to  the 
continuous  transform  of  (4.31)  and  (4.32).  For  the  above  DFT,  the  high¬ 
est  frequency  to  be  considered  Is  (N-l). Aw. 

The  fast  Fourier  transform  (FFT)  Is  an  efficient  method  for 
evaluating  DFT  based  on  the  numerical  algorithm  by  Cooley  and  Tukey 
(1965).  Basically,  the  FFT  makes  use  of  the  harmonic  nature  of  the  expo¬ 
nential.  By  setting  N  to  be  a  power  of  2  (l.e.,  2,  4,  8 . 1024, 

2048,.. etc.)  with  Integer  'm'  and  'n*  written  In  a  binary  form,  the  total 
computation  and  data  storage  Is  drastically  reduced.  There  are  many 
forms  of  FFT  available  which  differ  In  the  details  of  computer  coding, 
and  the  one  used  In  this  work  takes  advantage  of  a  complex  conjugate  sym¬ 
metry.  Thus,  for  the  real  Input  function  of  length  N,  the  routine  will 
return  only  half  of  the  transform  (from  zero  to  highest  frequency)  as 
complex  arrays  of  length  (N/2)+l  and  vice  versa  for  the  Inverse  FT. 

In  using  the  FFT,  the  values  of  basic  parameters  Involved  (num¬ 
ber  of  sampling  points  N,  time  Increment  At  and  period  Tp)  have  to  be 
selected  such  that  a  compromise  is  reached  between  the  accuracy  of 
results  and  the  cost  of  computations.  These  Investigations  will  be  pre¬ 
sented  In  Chapter  5. 
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In  summary,  the  frequency  domain  analysis  can  be  obtained  by: 
(1)  obtaining  the  Fourier  transform  of  the  loading  function,  (2)  multi¬ 
plying  this  Fourier  transform  by  the  transfer  function  of  the  system,  and 
(3)  obtaining  the  Inverse  Fourier  transform  of  the  product  to  get  the 
response  In  the  time  domain. 


4.3  NOTES  ON  COMPUTER  IMPLEMENTATION 


4.3.1  Integration 

To  form  matrices  G  and  H  In  the  formulation  In  Section  4.1,  the 
Integration  of  the  fundamental  solution  over  the  elements  as  defined  In 
(4.16)  needs  to  be  carried  out.  The  Integration  corresponding  to  the 
node  'k'  (where  the  virtual  load  Is  applied)  and  element  1  A*  has  no  prob¬ 
lem  as  long  as  k  t  t.  The  Integration  can  be  done  numerically  using  a 
four-point  standard  Gaussian  quadrature  Integration  (3x3  points  are  used 
In  the  3-D  case).  When  the  node  *  k'  belongs  to  element  *  £  *  (or  k  =  t), 
numerical  problems  arise  due  to  the  presence  of  a  singular  point  (at  r  = 
0)  In  the  Integral  domain.  In  this  case,  the  Integral  Hu  Is  zero  due  to 
the  orthogonality  of  the  normal  'n'  and  distant  'r'.  The  Integral  G^, 
however,  needs  different  treatment.  In  the  three-dimensional  case,  a 
series  expansion  of  the  exponential  function  was  used  and  the  Integral 
was  computed  analytically  (Dominguez,  1978).  In  the  two-dimensional 


Of 
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case,  the  numerical  Integration  was  still  used,  but  the  result  Is 
Improved  by  Introducing  a  correction. 

The  approach  to  compute  the  correction  starts  with  the  Investi¬ 
gation  of  the  function  U .,  as  shown  In  (4.3),  which  is  expressed  in  terms 

J  * 

of  modified  Bessel  functions  as  In  Equations  (4.7)  and  (4.8).  By  expand¬ 
ing  the  modified  Bessel  functions  Kn(Z)  Into  an  Infinite  series  (see 
Abramowltz  and  Stegun,  1965),  the  limit  forms  of  these  functions  are  tak¬ 
en  as  Z  approaches  zero  (for  this  case,  notes  that  Z  =  1«r/c,  where  c  can 

be  either  c,  ore  ) 
s  p' 

K0(Z)  =  -ln(§)  -  Y 

K,(Z)  =  \ 

(4.38) 

j  K-|(Z)  =  j2+  \  ln(^) 

k2(z)  =  h~  \ 

In  the  above  expressions,  %  is  an  Euler's  constant  which  has  no 
singularity  problem  In  numerical  Integration  and  can  be  taken  out.  Using 
the  relationship  In  (4.38),  the  limit  form  of  ♦  and  x  as  r  approaches 
zero  can  be  written  as 

*  =  (4'39) 
s  p 


x 


jt<VV2  - ” 


(4.40) 
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From  the  above  equation,  It  can  be  seen  that,  when  Is  Inte¬ 
grated  numerically,  the  only  term  that  contains  errors  due  to  singulari¬ 
ty  Is  Therefore,  only  Gjj  and  In  the  result  need  corrections. 
The  correction  which  has  to  be  added  to  G^  and  &££  Is  the  difference 
between  Integrating  (4.39)  analytically  and  numerically.  The  analytical 
Integration  of  (4.39)  can  be  written  as 


In  which  i  Is  the  length  of  the  boundary  element  under  consideration. 


4.3.2  Direction  of  Normal  and  Node  Arrangement 

The  problem  studied  Involves  formulations  for  both  cavities  and 
Inclusions,  which  are  the  Internal  and  external  surfaces  respectively. 
In  both  cases,  the  outer  normal  vector  n,  which  Is  necessary  In  defining 
the  direction  of  the  tractions,  Is  directed  outward  from  the  body.  The 
sequences  of  nodal  points  k,  l  In  each  element  need  to  be  defined  In  a 
manner  consistent  with  the  direction  of  the  outer  normal  for  that  partic¬ 
ular  element.  The  numbering  schemes  for  the  2-D  and  3-D  cases  used  In 
this  study  are  shown  In  Figure  4.4.  In  the  computer  program,  the  Input 
sequence  of  nodes  are  designed  to  follow  the  scheme  for  cavity.  The  com¬ 
puter  program  was  written  In  such  a  way  that  when  the  formulation  for 
Inclusions  is  needed,  It  will  reverse  the  nodal  sequence  temporarily  so 
that  the  outer  normal  n  will  be  In  the  opposite  direction. 


Two  Dimensional  Model 


Three  Dimensional  Model 


Figure  4.4  -  Node  Numbering  Scheme  for  External  and  Internal  Surfaces 
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4.4  DAMPING 

The  Internal  dissipation  of  energy  In  the  soil  is  believed  to  be 
due  to  Internal  friction  rather  than  to  viscoelastic  properties.  Con¬ 
stant  viscous  damping  produces  an  energy  loss  per  cycle  that  Increases 
linearly  with  frequency.  A  hysteretic  damping,  on  the  other  hand,  pro¬ 
duces  an  energy  loss  per  cycle  that  Is  frequency  independent,  but  depends 
on  the  amplitude  of  the  strains.  Experimental  evidence  (Hardin,  1965) 
indicates  that  the  second  case  Is  closer  to  the  true  behavior  of  soils. 
In  order  to  maintain  the  linearity  of  the  solution,  the  amplitude  depend¬ 
ence  Is  dropped,  using  what  Is  normally  called  a  linear  hysteretic  damp¬ 
ing.  In  the  frequency  domain  analysis,  this  type  of  damping  Is 
reproduced  by  assuming  a  complex  modulus  of  the  form  G(l+2iD),  where  D  Is 
the  damping  ratio.  The  terms  corresponding  to  Internal  soil  damping  are 
then  Included  In  the  complex  stiffness  matrix. 


CHAPTER  5 


PRELIMINARY  INVESTIGATIONS 

Since  finite  element  and  boundary  element  models  are  both  dis¬ 
crete,  the  results  are  expected  to  deviate  to  some  extent  from  the  the¬ 
oretical  solution  (If  It  Is  available)  depending  upon  the  nature  of  the 
model.  The  accuracy  of  the  results  can  be  Improved  by  making  the  dis¬ 
crete  model  become  as  close  as  possible  to  the  continuous  one.  This  can 
usually  be  done  by  choosing  a  very  small  value  of  the  discrete  parameters 
such  as  time  Increment,  size  of  element,  etc.  Of  course,  this  approach 
Is  not  always  practical  due  to  the  computational  costs.  Therefore,  It  Is 
necessary  to  Investigate  and  choose  a  suitable  value  of  these  parameters 
so  that  a  compromise  Is  reached  between  accuracy  and  cost  of  computation. 

In  this  chapter,  only  some  of  the  important  parameters  involved 
In  both  finite  element  and  boundary  element  formulations  are  Investi¬ 
gated  using  a  problem  similar  to  one  of  those  presented  In  Chapter  6. 
This  Is  a  case  of  free  field  wave  propagation  (no  Inclusions  Involved). 
The  target  Is  4  ft.  from  the  excitation  point  where  one  complete  cycle  of 
sinusoidal  force  Is  applied  In  a  direction  perpendicular  to  the  princi¬ 
pal  axis  of  wave  propagation,  l.e.,  only  a  shear  wave  Is  considered.  The 
frequency  of  excitation  f  is  100  cps.  The  medium  has  a  shear  wave 
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velocity  of  100  fps,  unit  weight  of  100  pcf,  a  Poisson's  ratio  of  0.25, 
and  no  material  damping. 

It  had  been  mentioned  In  Suddhlprakarn  (1983)  that  It  Is  simpler 
and  more  reliable  to  Interpret  the  results  If  the  motions  are  presented 
In  terms  of  a  displacement  rather  than  a  velocity  or  an  acceleration 
trace.  Therefore,  all  results  shown  In  this  work  will  be  In  the  form  of 
displacement  versus  time.  It  should  also  be  noted  that,  for  this  chap¬ 
ter,  the  results  are  still  presented  with  units  and  dimensions.  In  the 
next  chapters,  the  dimensionless  form  will  be  used. 


5.1  FINITE  ELEMENT 


Because  the  appropriate  value  of  the  time  step  of  Integration 
used  In  the  finite  element  model  can  be  easily  determined  using  the  meth¬ 
ods  outlined  In  Section  3.4,  they  will  not  be  discussed  here.  The  only 
preliminary  Investigation  needed  for  the  finite  element  method  In  this 
work  Is  to  determine  the  size  of  the  domain  used  In  the  model  In  order  to 
represent  properly  an  Infinite  medium. 

Figure  5.1  shows  the  finite  element  model  with  the  size 
expressed  relative  to  the  distance  'd1  between  the  point  of  excitation 
'A*  and  the  target  'B'.  In  this  case,  the  boundaries  may  have  effects  on 
the  motion  at  both  the  excitation  point  and  the  target.  There  are  three 
reflected  waves  Involved  In  this  case:  the  first  one  Is  a  shear  wave 
reflected  from  point  'O'  at  the  left  boundary  passing  'A'  to  the  target 
' B* ,  the  second  one  Is  also  a  shear  wave  reflected  from  the  right  bounda- 


Soil  Properties:  v  *  0.25 

Y  =  100  pcf 
cs  =  100  fps 

Frequency  of  Excitation  =  100  Hz 

Figure  5.1  -  Finite  Element  Model  Representing  an  Infinite  Medium 


ry  at  '£'  back  to  the  target  '8',  and  the  third  one  Is  a  compression  wave 
that  takes  the  shortest  route  to  the  bottom  boundary  at  point  'C'  and 
reflects  back  to  the  target  point  'B1 . 

Although  the  restraining  effect  and  the  reflected  wave  effects 
both  happen  together,  It  Is  possible  to  differentiate  between  them 
because  the  time  Involved  In  the  travelling  of  reflected  waves  can  be 
estimated.  Once  the  arrival  time  of  reflected  waves  Is  put  beyond  the 
time  span  of  our  Interest,  the  difference  In  the  results  If  any  exists, 
will  be  caused  solely  by  the  restraining  effects  from  the  boundaries. 

In  this  section,  the  effects  from  each  boundary  are  Investigated 
separately.  By  keeping  the  other  two  boundaries  far  away  to  ensure  that 
they  cause  no  effects,  the  boundary  under  consideration  can  be  located  In 
such  a  way  that  there  Is  no  further  significant  changes  In  the  result 
when  this  boundary  Is  moved  away  from  the  travelling  path  of  the  direct 
wave. 

The  effect  from  the  left,  right  and  bottom  boundaries  are  shown 
In  Figures  5.2,  5.3,  and  5.4  respectively.  In  each  figure,  there  are 
four  separate  plots  of  displacement  versus  time  for  different  locations 
of  the  boundary  under  study.  All  notations  related  to  the  size  of  the 
domain  (Ml,  N2,  and  N3  in  this  case)  are  defined  in  Figure  5.1.  The  the¬ 
oretical  arrival  time  for  the  plane  shear  wave,  the  compression  wave,  the 
reflected  shear  wave,  and  the  reflected  compression  wave  are  denoted  by 
the  vertical  lines  marked  S,  P,  SR,  and  PR  respectively.  In  all  results, 
there  Is  an  early  motion  occurring  at  the  time  of  arrival  of  the  com1 
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presslon  wave  which  Is  due  to  the  nature  of  the  point  force  excitation 
(see  Suddhlprakarn,  1983). 

In  Figure  5.2,  the  effects  of  the  left  boundary  are  Investigated 
by  setlng  N2  and  N3  constant  at  1.75.  The  results  for  the  case  N1  of 
0.25,  0.50,  0.75,  and  1.0  are  shown  In  Figures  5.2a,  5.2b,  5.2c,  and  5. 2d 
respectively.  The  first  thing  to  note  Is  that  the  shapes  of  the  dis¬ 
placement  curves  before  the  arrival  of  the  reflected  wave  are  the  same 
for  all  four  cases.  The  only  difference  Is  that  the  amplitude  Is  slight¬ 
ly  larger  in  the  case  when  the  left  boundary  Is  closer  to  the  point  of 
excitation,  especially  In  the  case  N1  =  0.25.  For  N1  greater  than  0.50, 
the  amplitude  remains  constant. 

In  Figure  5.2a,  the  arrival  of  the  shear  wave  reflected  from  the 
left  boundary  Is  evident  by  a  starting  downward  motion  of  quite  large 
amplitude.  This  Is  In  contrast  to  the  case  of  the  direct  wave  which 
starts  with  an  upward  motion  conforming  to  the  excitation  force.  The 
opposite  nature  of  the  direct  and  reflected  shear  waves  Is  reasonable  If 
one  considers  the  zero  displacements  at  the  fixed  boundaries. 

For  N1  of  0.50  and  higher  as  in  Figures  5.2b,  5.2c,  and  5. 2d,  the 
arrival  time  of  the  reflected  shear  wave  Is  outside  the  time  span  pre¬ 
sented  and  only  the  restraining  effect  from  the  boundary  shows  up.  It  Is 
readily  seen  that,  without  the  interference  of  the  reflected  wave,  the 
motions  of  all  three  cases  In  Figures  5.2b  to  5. 2d  Is  practically  the 
same  except  that  Figure  5.2b  shows  a  slightly  larger  amplitude  of  the 
negative  peak.  In  the  cases  of  N1  =  0.75  and  1.00,  the  displacements  are 
Identical  except  for  the  decaying  part  of  the  motion  after  t  =  0.06  sec. 
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which  oscillates  more  uniformly  In  the  case  of  N1  =  1.00.  This  motion 
could  continue  for  a  long  time  especially  In  a  system  without  material 
damping  like  this  one;  therefore,  It  should  not  be  a  major  factor  In 
choosing  N1  to  be  used  as  long  as  the  significant  part  of  the  motion  Is 
accurate  enough.  For  this  case,  the  value  of  N1  =  1.00  seems  to  give  a 
satisfactory  result. 

Using  N1  =  1.00  and  keeping  N3  constant  at  1.75,  the  effects  of 
the  right  boundary  can  be  studied  In  the  same  manner.  The  results  for 
the  cases  N2  of  0.25,  0.50,  1.00,  and  1.50  are  shown  In  Figures  5.3a, 
5.3b,  5.3c,  and  5.3d  respectively.  The  results  from  Figures  5.3a  and 
5.2a  are  almost  Identical,  which  shows  that  the  shear  waves  reflected 
from  the  left  and  the  right  boundary  have  the  same  effects.  By  comparing 
Figure  5.3  with  Figure  5.2,  It  Is  apparent  that  the  effects  of  the  right 
boundary  show  the  same  pattern  discussed  In  the  case  of  the  left  bounda¬ 
ry.  Figure  5.3c  and  Figure  5.3d,  for  N2  of  1.0  and  1.5  respectively, 
show  very  similar  results  except  for  the  last  cycle  of  the  decaying  por¬ 
tion  which  Is  considered  unimportant.  The  further  Increase  of  N2  from 
1.50  (as  In  Figure  5.3d)  to  1.75  (as  In  Figure  5. 2d)  does  not  Improve  the 
results.  From  these  results,  It  appears  that  the  suitable  value  for  N2 
would  range  from  1.0  to  1.75. 

In  Figure  5.4,  the  effects  of  the  lower  boundary  are  Investi¬ 
gated.  This  time  the  left  and  right  boundaries  are  far  away  (N1  =  1.0, 
N 2  =  1.75).  Figures  5.4a,  5.4b,  5.4c,  and  5.4d  show  the  results  for  the 
cases  of  N3  equals  to  0.25,  0.50,  1.0,  and  1.50  respectively.  In  all 
four  cases,  the  reflected  compression  waves  arrive  at  the  target  within 
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the  time  span  presented.  This  compression  wave  Is  also  reflected  from 
the  top  boundary  which  Is  not  shown  In  Figure  5.1  since  the  model  takes 
advantage  of  symmetry  and  antisymmetry  as  mentioned  before.  In  Figure 
5.4a  where  N3  Is  0.25  (which  means  that  the  dimension  N3.d  In  Figure  5.1 
Is  the  same  as  the  wavelength),  the  motion  Is  extremely  complicated  by 
the  reflected  compression  wave  which  arrives  even  earlier  than  the 
direct  shear  wave.  In  Figure  5.4b  where  N3  Is  Increased  to  0.50,  the 
reflected  compression  wave  Is  still  ahead  of  the  direct  shear  wave  but 
the  motion  In  this  case  Is  slightly  better.  For  N3  of  1.0  and  1.50  as 
shown  In  Figures  5.4c  and  5.4d  respectively,  the  reflected  compression 
wave  arrives  after  the  direct  shear  wave  and  the  motions  are  more  accu- 
rate  except  after  the  arrival  of  the  reflected  compression  wave  which 
still  affects  the  motion.  The  comparison  of  Figure  5.4d  and  Figure  5. 2d, 
with  N1  equals  to  1.50  and  1.75  respectively,  show  that  the  motion  before 
the  arrival  of  the  reflected  compression  wave  Is  identical,  and  the  same 
conclusion  can  also  be  drawn  from  Figures  5.4c  and  5.4d.  This  Indicates 
that  the  major  factor  to  be  considered  on  selecting  the  value  of  N3  is  to 
keep  the  reflected  compression  waves  beyond  the  time  span  of  Interest  and 
for  this  case,  a  value  of  N3  of  1.75  would  be  appropriate. 

In  summary,  to  represent  an  Infinite  domain  by  the  finite  ele¬ 
ment  model,  as  shown  In  Figure  5.1,  the  main  consideration  Is  to  keep  the 
reflected  waves  out  of  the  time  span  under  consideration.  This  Is  usual¬ 
ly  good  enough  for  a  bottom  boundary  (and  a  top  boundary  as  well).  For 
the  right  and  left  boundaries,  they  still  need  to  be  extended  slightly 
farther  so  that  the  motion  in  the  decaying  zone  remains  practically  the 
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same.  For  the  case  In  this  chapter,  the  value  of  Nl,  N2,  and  N3  of  1.0, 
1.50,  and  1.75,  respectively,  seem  to  be  a  good  compromise  between  accu¬ 
rate  results  and  minimum  cost  of  computation. 


5.2  BOUNDARY  ELEMENT 

The  preliminary  studies  on  the  boundary  element  model  presented 
in  this  chapter  deal  almost  entirely  with  the  parameters  used  for  the 
fast  Fourier  transform  except  for  the  first  subsection  which  discusses 
the  methods  to  handle  the  problem  of  the  transfer  function  at  zero  fre¬ 
quency. 


5.2.1  Transfer  Function  at  Zero  Frequency 

As  mentioned  In  Chapter  4,  In  a  frequency  domain  analysis,  a  sol¬ 
ution  or  a  transfer  function  has  to  be  solved  for  different  frequencies 
of  a  unit  point  excitation  from  zero  to  the  highest  frequency  required  In 
the  fast  Fourier  transform.  A  problem  arises  when  attempting  to  solve 
the  displacements  from  the  formulation  shown  In  Equations  (4.25)  and 
(4.26)  at  zero  frequency,  because  the  term  *  and  x  are  singular  at  u  =  0 
as  can  be  seen  from  Equations  (4.5)  to  (4.8).  This  causes  the  Green's 
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function  Ujj  and  T^  to  be  discontinuous  at  zero  frequency. 

Figure  5.5  Is  the  typical  shape  of  the  transfer  function  versus 
frequency  curve  In  our  problem.  The  transfer  functions  plotted  on  the 


Transfer  Function 


y-axls  are  the  absolute  values  since  they  are  complex  quantities.  At 
large  frequencies,  the  curve  Is  almost  horizontal,  then  it  slopes  up  gra¬ 
dually  when  the  frequency  keeps  decreasing.  At  small  frequencies,  the 
tranfer  function  Increases  rapidly  and  approaches  infinity  as  the  fre¬ 
quency  approaches  zero  as  shown  in  the  figure. 

Three  methods  were  tried  to  handle  this  problem,  and  they  are 
shown  in  Figure  5.6.  The  first  case  shown  in  Figure  5.6a  is  using  the 
static  solution  at  zero  frequency  Instead  of  the  dynamic  solution  as  in 
Equations  (4.3)  and  (4.4). 

Since  the  dynamic  solution  does  not  approach  the  static  one  when 
the  frequency  approaches  zero,  the  curve  of  the  transfer  function  near 
zero  frequency  can  have  different  shape  depending  upon  the  size  of  the 
frequency  increment  Af.  Two  possible  shapes  are  shown  in  Figure  5.6a  as 
either  'ade'  for  large  Af,  or  'abcde'  for  small  Af. 

In  the  second  case,  the  static  solution  is  used  as  an  upper  bound 
as  shown  in  Figure  5.6b.  In  this  case,  the  shape  of  the  transfer  func¬ 
tion  curve  Is  independent  of  Af.  In  the  third  case,  the  dynamic  solution 
is  used  throughout,  but  at  zero  frequency,  an  assumed  value  that  is  small 
relative  to  Af  Is  used  Instead  of  zero  to  avoid  numerical  problem.  The 
transfer  function  in  this  case  Is  also  Independent  of  Af  except  that  the 
starting  point  of  the  curve  will  depend  on  the  value  assumed  for  zero 
frequency. 

The  above  three  methods  are  compared  by  analysing  the  free  field 
wave  propagation  problem.  For  the  range  of  the  Fast  Fourier  transform 
parameters  that  are  going  to  be  used  in  this  work,  the  solutions  from  all 


three  approaches  are  identical  and  will  not  be  presented  here.  The  rea¬ 
son  for  these  same  results  Is  that  the  transfer  functions  In  all  three 
cases  only  differ  at  a  few  points  near  a  zero  frequency.  These  few 
points,  in  spite  of  their  large  values,  contribute  very  little  to  the 
Fourier  transform. 

If  the  first  or  the  second  method  above  is  going  to  be  used  and 
the  problems  other  than  a  free  field  case  need  to  be  analysed,  the  stat¬ 
ic  solution  has  to  be  Incorporated  consistently  into  the  different  parts 
of  the  computer  program,  such  as  the  calculation  of  the  matrices  G  and  H. 
This  involves  extra  works  which  is  unnecessary  since  the  third  approach 
seems  to  produce  equally  good  results  without  any  significant  changes  in 
the  computer  program.  Therefore,  the  third  approach  will  be  used 
throughout  this  work. 


5.2.2  Fast  Fourier  Transform  Parameters 

The  basic  concept  of  the  Fourier  transform  and  a  frequency 
domain  analysis  were  discussed  In  Section  4.2.  The  analysis  In  the  fre¬ 
quency  domain  Involves  three  basic  steps:  (1)  Obtaining  the  Fourier 
transform  of  the  loading  function,  (2)  Multiplying  this  Fourier  trans¬ 
form  by  the  transfer  function  of  the  system,  and  (3)  Obtaining  the 
inverse  Fourier  transform  of  the  results  of  step  (2)  to  get  the  time 
response.  The  parameters  Involved  In  the  fast  Fourier  transform  (FFT) 
will  be  discussed  by  following  these  three  logical  steps. 


•  V-. 
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In  step  (1),  to  obtain  the  FFT  of  a  nonperiodic  loading  function, 
the  period  of  this  loading  function,  which  Is  considered  to  be  infinite 

In  a  classical  Fourier  transform,  has  to  be  assumed  as  a  finite  value  T  . 

P 

The  loading  function  used  In  this  study  Is  shown  In  Figure  5.7a.  It  Is 
apparent  that  the  loading  function  shown  will  represent  exactly  a  single 
cycle  sinusoidal  force  If  the  period  Tp  Is  extended  to  infinity.  This 
Is,  of  course.  Impossible  numerically.  In  practice,  however,  it  is 
acceptable  to  select  a  period  Tp  sufficiently  large  so  that  the  effects 
of  the  spurious  repetitive  loading  (shown  as  the  dotted  curve  in  Figure 
5.7a)  are  negligible. 

In  FFT,  this  continuous  loading  function  is,  then,  sampled  at 
discrete  time  interval  At  throughout  the  period  Tp  as  shown  In  Figure 
5.7a.  The  total  number  1 N 1  of  the  sampling  points  and  the  time  increment 
At  will  decide  the  value  of  the  period  T  .  Furthermore,  they  will  also 
determine  the  frequency  Increment  Af,  the  maximum  frequency  fm  ,  and 
the  number  of  frequency  points  at  which  the  transfer  function  has  to 
be  calculated.  The  relationship  of  these  parameters  can  be  written  as 


Tp  =  NAt 
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It  should  be  noted  that  the  number  of  sampling  points  N  has  to  be 
a  power  of  '2'  as  mentioned  in  Section  4.2  for  the  FFT  algorithm. 

It  can  be  seen  from  the  expressions  above  that  the  fundamental 
parameters  involved  are  only  N  and  At  from  which  the  other  parameters  can 
be  obtained.  Therefore,  Investigations  will  be  made  first  for  N  and  At. 

(a)  Effects  of  N  Different  values  of  N  shown  In  Figure  5.7b  are  used  with 

the  value  of  At  fixed  at  0.0005  seconds.  The  period  T  for  each  case 

P 

relative  to  the  period  T$  of  the  sinusoidal  part  of  the  loading  function 
are  also  tabulated.  The  results  for  each  value  of  N  are  shown  In  Figures 
5.8a  to  5. 8d.  In  Figures  5.8a  and  5.8b  where  N  Is  1024  and  512  respec¬ 
tively,  there  is  no  difference  In  the  results.  When  N  is  reduced  to  256 
as  in  Figure  5.8c  (where  Tp/Ts  Is  now  12.8),  the  positive  and  negative 
peaks  are  still  the  same  as  In  the  first  two  figures  but  the  beginning 
part  of  the  curve  starts  to  show  a  very  small  negative  value  Instead  of 
zero.  It  Indicates  that  the  period  Tp  is  not  large  enough,  and  conse¬ 
quently,  the  effects  of  the  repetitive  loading  (shown  as  a  dotted  curve 
on  the  left  in  Figure  5.7a)  become  significant.  The  result,  however,  Is 
still  acceptable  unlike  the  case  of  N  =  128  shown  In  Figure  5.8d  In  which 
the  beginning  part  of  the  curve  shows  larger  deviations  from  zero  (even 
though  the  other  part  of  t..«  result,  still,  seems  to  be  satisfactory). 

Figure  5.8d  also  shows  that  the  curve  terminates  before  the  end 
of  the  time  scale;  therefore,  the  time  span  under  study  has  to  be  consid¬ 
ered  as  well  In  choosing  the  value  of  N  and  At. 


All  In  all,  when  the  number  of  the  sampling  points  N  decreases 
such  that  the  period  Tp  Is  closer  to  the  period  T$  of  the  sinusoidal  part 
of  the  loading  function,  the  effects  of  the  repetitive  loading  start  to 
show  up.  In  other  words,  the  quiet  duration  between  point  T$  and  Tp  In 
Figure  5.7a  has  to  be  long  enough  so  that  the  transient  part  from  the 
previous  loading  cycle  completely  disappears. 

Studies  were  also  conducted  for  the  case  of  closer  distances 
between  the  excitation  point  and  the  target,  and  different  period  T  . 
The  results,  though  not  presented  here,  show  that  acceptable  accuracy  Is 
obtained  If  the  ratio  Tp/Ts  Is  not  lower  than  a  value  of  about  10. 

(b)  Effect  of  At  Different  values  of  At  are  Investigated  with  the  value 
of  N  fixed  at  2048.  The  only  two  cases  presented  are  tabulated  In  Figure 
5.7c  together  with  the  corresponding  values  of  T  .  In  both  cases,  the 
period  ^  Is  long  enough  so  that  no  transient  from  the  previous  cycle 
appears.  The  results  for  the  case  At  of  0.00025  seconds  and  0.001  sec¬ 
onds  are  shown  In  Figures  5.9a  and  5.9b  respectively.  The  general  shape 
of  the  curve  Is  pretty  much  the  same  In  both  cases  except  that  for  the 
larger  At,  the  resolution  of  the  curve  may  not  be  sufficient  to  reproduce 
some  sharp  peaks  as  shown  In  Figure  5.9b. 

Up  to  this  point,  the  only  conclusions  that  can  be  drawn  about 
the  selection  of  At  Is  that  It  should  be  sufficiently  small  In  order  to 
get  smooth  results.  In  addition,  a  small  value  of  At  Is  also  necessary 
to  reproduce  accurately  the  loading  function.  On  the  other  hand,  At 


Figure  5.9  -  Effects  of  the  Time  Increment  at 
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should  not  be  so  small  that  It  reduces  the  value  of  the  period  causing 
the  transient  to  show  up  again. 

In  step  (1),  the  N  points  sampled  from  the  loading  function  are 

real  quantities  in  the  time  domain.  After  taking  a  Fourier  transform, 

the  results  consist  of  N^  complex  quantities  In  the  frequency  domain, 

from  zero  to  the  highest  frequency  f  ,  with  an  increment  of  Af.  In 

max 

step  (2),  each  of  these  complex  points  Is  then  multiplied  by  the  values 
of  the  transfer  function  determined  at  each  corresponding  frequency. 

Since  the  calculation  of  the  transfer  function  Is  time  consuming 
especially  for  the  case  with  inclusions,  It  Is  usually  not  possible  to 
compute  the  transfer  function  at  every  required  frequency  If  the  compu¬ 
tation  cost  has  to  be  kept  within  a  practical  limit.  Normally,  only  some 
of  the  values  of  the  transfer  function  need  to  be  calculated  while  the 
remaining  ones  can  be  defined  In  different  way  depending  on  the  procedure 
used.  In  this  study,  two  procedures:  truncation,  and  interpolation  are 
investigated.  Both  procedures  Involve  the  parameters  f  ,  Af,  and 
shown  In  Equation  (5.1). 

(c)  Effects  of  truncation  Figure  5.10a  shows  typical  transfer  functions 
at  each  frequency  Increment  Af,  from  zero  to  the  highest  frequency  f  . 

max 

In  the  truncation  procedure,  the  transfer  functions  are  calculated  at 
each  frequency  Increment  Af  up  to  the  truncated  frequency  f^  ,  beyond 
which  they  are  assumed  to  be  constant  for  case  1,  and  are  defined  as  zero 
for  case  2,  as  shown  In  Figures  5.10b  and  5.10c  respectively. 


Transfer  Function  Transfer  Function  Transfer  Function  Transfer  Function 


(a)  All  Points  Calculated 


(b)  Truncation,  Case  1 


Frequency 

(c)  Truncation,  Case  2 


Figure  5.10  -  Truncations  and  Interpolation  of  the 
Transfer  Function 


Using  the  truncation  procedure  case  1,  the  results  for  different 
value  of  ftr  are  shown  In  Figures  5.11a  to  5. lid.  In  each  case,  the  rel¬ 
ative  value  of  ftr  with  respect  to  fmax  and  f  (the  frequency  of  the 
sinusoidal  part  of  the  loading  function)  are  also  shown.  In  Figure 
5.10a,  the  transfer  functions  are  calculated  at  every  point  without  any 
truncation,  resulting  In  a  very  smooth  curve.  In  Figure  5.11b  where  f^ 
equals  l/2f  (5f  ),  there  are  some  noticeable  vibrations  at  the  begin- 
nlng  of  the  curve  and  the  overall  result  Is  slightly  less  smooth  than  In 
Figure  5.11a;  however,  the  result  Is  still  acceptable.  In  Figure  5.11c, 
where  ftr  equals  l/4fmax  (and  2.5fs),  many  Irregularities  are  apparent 
although  the  peak  value  Is  still  accurate.  In  Figure  5. lid,  the  results 
for  two  different  values  of  ftr  are  plotted  together,  the  solid  curve 

corresponding  to  f.„  =  l/8f  „  and  the  dotted  curve  to  l/16f  .  In  both 

tr  max  max 

cases  the  curves  are  extremely  Irregular  and  the  results  lack  sufficient 
accuracy. 

The  corresponding  set  of  results  using  the  truncation  procedure 

case  2  are  shown  In  Figure  5.12.  As  before,  Figure  5.12a  shows  the 

result  without  any  truncation  just  for  comparison.  By  comparing  the 

results  between  the  truncation  case  1  and  case  2  for  each  corresponding 

value  of  ftr.  It  Is  apparent  that  the  results  from  both  cases  are  the 

same  except  the  beginning  part  of  the  curve  which  Is  now  flatter  In  the 

truncation  case  2,  even  when  r  Is  as  low  as  l/4f  „  (Figure  5.12c).  In 

tr  max 

Figure  5.12d,  even  though  the  high  frequency  components  are  not  enough  to 
produce  accurate  results,  the  curves  still  show  more  uniform  motions  at 
the  beginning  part  comparing  with  those  of  Figure  5. lid. 
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The  results  of  these  studies  Indicate  that  to  obtain  a  reason¬ 
able  accuracy,  making  zero  the  transfer  functions  after  f  would  be  a 
better  choice  than  assuming  constant  values.  When  the  calculated  trans¬ 
fer  functions  are  used,  the  results  (Figure  5.12)  Indicate  that  the 
transfer  functions  at  higher  frequencies  have  small  effects  on  the 
results.  The  frequency  beyond  which  transfer  functions  become  less 
effective  is  believed  to  be  related  to  the  frequency  of  the  sinusoidal 
excitation  f  ,  and  this  was  confirmed  by  further  studies  not  shown  here 
using  different  values  of  f$  and  f.  .  All  these  studies  show  that  the 
result  will  be  accurate  enough  as  long  as  f^  Is  more  than  five  times  f  . 

(d)  Effects  of  Interpolation  The  Interpolation  procedure,  represented 
graphically  In  Figure  5.10d,  Involves  the  calculation  of  the  transfer 
functions  from  zero  to  the  highest  frequency,  but  at  a  frequency  Incre¬ 
ment  different  from  Af  defined  In  Equation  (5.1).  The  frequency  Incre¬ 
ment  used  In  the  Interpolation  method,  denoted  as  'DFR1,  Is  arbitrarily 
chosen,  larger  than  Af,  and  not  necessarily  a  multiple  of  Af ;  therefore, 
the  number  of  frequencies  at  which  the  transfer  functions  are  calculated 
(NFR)  Is  less  than  N^  and  these  frequencies  may  not  coincide  with  the 
frequencies  where  the  transfer  functions  are  required.  The  transfer 
function  at  each  required  frequency  Is  determined  using  a  linear 
Interpolation  as  shown  In  Figure  5 . lOd . 

Different  combinations  of  the  number  of  frequencies  NFR  and  the 
frequency  Increment  DFR  were  Investigated.  All  combinations  were  chosen 
such  as  to  cover  the  complete  range  of  frequencies  from  zero  to  fMW,  and 
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the  results  are  shown  in  Figures  5.13a  to  5.13d.  In  Figure  5.13a  when 
NFR  in  1001  and  DFR  is  1,  the  results  are  almost  the  same  as  in  the  case 
of  Figure  5.12a  where  no  truncation  is  used;  therefore,  both  results  are 
smooth  with  a  slightly  lower  peak  in  the  case  of  Figure  5.13a.  When  DFR 
increases  to  2,  and  NFR  Is  reduced  to  501  as  In  Figure  5.13b,  the  shape 
of  the  curve  is  the  same  but  the  peak  is  lower.  This  decrease  in  the  peak 
amplitude  Is  more  pronounced  when  less  Interpolation  points  are  used  as 
evidenced  in  Figures  5.13c  and  5.13d.  In  the  latter  figure,  the  shape  of 
the  curves  is  also  badly  deformed. 

Although  the  shapes  of  the  displacement  curves  in  the  case  of 
interpolation  are  almost  the  same  throughout  (except  when  DFR  is  large), 
the  amplitude  of  the  peak  Is  reduced  In  comparison  to  the  results  from 
the  truncation  procedure  that  use  a  comparable  number  of  points  for  the 
calculation  of  the  transfer  function  (see  Figures  5.13b  and  5.12b  for 
example).  It  appears  .therefore,  that  truncation  of  the  transfer  func¬ 
tion  provides,  for  the  same  computational  effort,  a  better  solution  than 
interpolation. 


5.2.3  Summary  on  Parameters  of  the  Frequency  Domain  Analysis 

For  the  frequency  domain  analysis  to  be  useful,  the  mumber  of 
points  at  which  the  transfer  function  Is  calculated  should  be  within  the 
range  of  100  to  200  In  order  to  maintain  the  cost  of  computation  as  low 
as  possible.  The  method  of  truncation  case  2  seems  to  be  the  best  method 
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to  keep  the  number  of  points  within  this  range  without  sacrificing  the 

accuracy  of  the  results. 

From  the  above  considerations,  the  following  conclusion  can  be 

drawn 

1.  At  zero  frequency,  where  a  singularity  problem  occurs,  It  Is  easier 
and  more  consistent  to  compute  the  transfer  function  at  a  very  small 
frequency  Instead  of  using  the  correct  static  values. 

2.  The  truncation  method  case  2  proves  to  be  the  most  efficient  way  to 
reduce  the  number  of  the  transfer  functions  that  need  to  be  calcu¬ 
lated. 

3.  The  time  Increment  At  should  be  small  enough  to  give  smooth  results 
with  well-defined  peaks,  and  to  reproduce  accurately  the  loading 
function. 

4.  The  combination  of  At  and  the  number  of  sampling  points  N  should 
result  In  a  period  Tp  that  Is  long  enough  to  provide  sufficient  quiet 
duration  beyond  the  sinusoidal  part  of  the  loading  function  so  that 
the  effects  of  the  repetitive  loading  disappear.  The  period  Tp 
should  be  more  than  ten  times  T$  (the  period  of  the  sinusoidal  part 
of  the  loading  function).  In  addition,  the  combination  of  At  and  N 
should  produce  a  total  time  long  enough  to  cover  the  time  span  of 
Interest. 
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5.  The  truncated  frequency  f^r  should  be  selected  such  that  there  are 
enough  high  frequency  components  In  the  Fourier  transform  to 
reproduce  the  wave  frequency.  For  this  study,  ft  should  be  more 
than  five  times  f  the  frequency  of  the  sinusoidal  part  of  the  load¬ 
ing  function. 

These  preliminary  studies  were  Intended  to  determine  the  effect 
of  various  parameters  on  the  accuracy  of  the  frequency  domain  solution 
using  the  FFT.  Two  additional  points  that  must  be  considered  are  the 
distance  from  the  excitation  to  the  target,  and  the  size  of  the  boundary 
elements.  Clearly  If  one  wants  to  determine  the  time  of  arrival  of  the 
first  wave  at  the  target  with  an  error  less  than  n  percent,  the  time  step 
At  should  not  be  larger  than  (n/100).d/c  where  d  Is  the  distance,  and  c 
the  wave  propagation  velocity.  To  reproduce  properly  the  transfer  func¬ 
tion  at  a  frequency  f,  the  size  of  the  boundary  elements,  similarly  to 
the  part  made  for  the  finite  element  mesh,  should  not  be  larger  than  X/4 
to  X/8  where  X  Is  the  wavelength  c/f. 
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CHAPTER  6 

EFFECT  OF  AN  INCLUSION  ON  ATTENUATION 

The  studies  of  wave  attenuation  were  performed  by  observing  the 
peak  amplitude  and  the  shape  of  the  displacement  versus  time  curve  at 
various  distances  from  the  point  of  excitation.  Analyses  were  made  (both 
In  two  and  three  dimensions)  for  the  case  with  Inclusions.  For  the  boun¬ 
dary  element  model  developed  In  this  work,  all  four  physical  properties 
of  soil  and  Inclusions— which  are  shear  wave  velocity,  mass  density, 
Poisson's  ratio,  and  damping  ratio — can  be  varied  In  any  combination  as 
required.  By  approximating  the  curved  surface  with  several  straight 
elements,  different  shapes  and  sizes  of  Inclusions  can  be  modeled  with 
sufficient  accuracy.  In  addition,  the  cluster  of  Inclusions  can  be 
arranged  In  any  configuration.  The  only  limitation  Is  In  the  total  num¬ 
ber  of  boundary  elements  due  to  the  size  of  the  resulting  matrix  that  can 
be  handled  In  core  by  the  computer. 

The  possible  combination  of  the  parameters  Involved  Is  so  large 
that  a  complete  parametric  study  Is  virtually  Impossible.  Only  a  few 
typical  cases  for  the  shear  wave  will,  therefore,  be  studied  In  this 
work.  In  this  chapter,  results  of  comparative  studies  between  the  free 
field  case  and  the  case  with  a  square  inclusion  (a  cube  for  the 


three-dimensional  case)  having  different  sizes  In  relation  to  the  wave¬ 
length  X.  In  each  case,  different  values  of  soil  damping  are  used. 

As  presented  In  Stokoe  et  al  (1980),  the  finite  element  method 
(FEM)  was  Implemented  for  the  wave  propagation  problem  In  the  trl axial 
cube  containing  soil  and  rigid  Inclusions,  The  method  Is  more  appropri¬ 
ate  because  of  the  finite  domain,  while  the  boundary  element  method 
(BEM),  on  account  of  a  need  to  discretize  the  outer  surface  of  the  cube 
as  well,  Is  very  expensive  for  that  particular  problem.  On  the  other 
hand,  for  an  Infinite  medium  as  In  this  work,  the  FEM  Is  not  as  economic 
and  flexible  as  the  BEM;  therefore,  It  Is  used  only  to  compare  the 
results  It  provides  with  those  of  the  boundary  element  solution.  The 
comparison  of  results  Is  presented  In  the  first  section  of  this  chapter 
with  all  results  Involving  the  attenuation  of  the  shear  wave  In  the  free 
field. 

All  results  are  presented  In  dimensionless  form.  The  notation, 
unless  defined  where  It  Is  used,  will  be  shown  In  the  list  of  symbols. 
The  most  frequently  used  forms  are  a  dimensionless  displacement  uG/F  and 
a  dimensionless  time  c$t/d,  where  u  Is  the  displacement  at  the  target,  G 
Is  the  shear  modulus  of  the  soil,  F  Is  the  peak  amplitude  of  the  excita¬ 
tion  (force  per  unit  length  for  the  Z-D  case),  cs  is  the  shear  wave 
velocity  of  the  soli,  t  Is  the  elapsed  time,  and  d  Is  the  distance 
between  the  source  and  the  target.  In  this  way,  the  theoretical  time  of 
arrival  of  the  shear  wave  In  the  free  field  corresponds  to  a  value  of 


6.1  FREE  FIELD  ATTENUATION 


Due  to  the  limited  capability  of  the  finite  element  model,  com¬ 
parison  of  the  results  from  the  FEW  and  the  BEM  can  be  made  only  for  the 
two  dimensional  case  without  material  damping  as  will  be  presented  In  the 
following  subsection. 


6.1.1  Two-Dimensional  Casa  with  tha  FEM  Results 

The  displacement  versus  time  at  different  distance  ratio  d/X, 
using  the  FEM  and  the  BEM  are  shown  In  Figures  6.1  and  6.2  respectively. 
For  the  finite  element  result  shown  In  Figure  6.1,  the  size  of  the  domain 
used  cannot  be  further  enlarged  to  accommodate  the  whole  range  of  d/X 
without  a  substantial  Increase  In  the  computation  cost.  Consequently, 
the  reflected  waves  will  be  allowed  to  appear  In  the  time  span  presented 
(after  the  completion  of  the  direct  wave)  provided  that  they  do  not 
affect  the  arrival  time  and  peak  amplitude.  Only  the  case  of  d/X  =  10, 
where  both  the  target  and  the  point  of  excitation  are  too  close  to  the 
boundary,  the  result  seems  to  be  unreliable  and  will  not  be  used. 

The  results  from  both  FEM  and  BEM  show  the  small  early  motion 
starting  at  the  arrival  time  of  the  plane  p-wave.  The  peak  amplitude  of 
these  early  motions  attenuates  very  fast,  as  can  be  seen  from  both  fig¬ 
ures.  At  a  large  value  of  d/X  when  the  cylindrical  wavefront  can  be 
approximated  as  a  plane,  these  early  motions  are  barely  discernible. 
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The  difference  between  the  results  from  both  methods  is  clearly 
seen  from  Figures  6.1  and  6.2.  The  first  point  to  note  is  that  the 
amplitude  of  the  positive  peak,  from  the  FEM  is  lower  and  attenuates  at  a 
faster  rate  than  that  from  the  6EM.  Both  positive  and  negative  peaks  in 
the  FEM  results  are  almost  the  same  in  magnitude  and  occur  slightly  after 
the  corresponding  peaks  In  the  BEM  results  whose  negative  peak  is  approx¬ 
imately  one  half  of  the  positive  one.  Secondly,  the  motion  after  the 
loading  cycle  In  the  BEM  decreases  asymptotically  while  In  the  FEM,  the 
motion  is  oscillatory  with  a  decay  In  amplitude  with  time  similar  to  what 
could  be  expected  from  a  lightly  damped  single  degree  of  freedom  system. 
Finally,  the  overall  shape  of  the  displacement  versus  time  curve  from  the 
FEM  is  smoother.  This  Is  evident  especially  at  the  point  of  the  first 
arrival  of  the  S-wave  where  the  displacement  from  the  FEM  solution 
Increases  less  abruptly. 

It  should  be  noted  that  the  two-point  excitation,  simulating  the 
condition  in  the  finite  element  model.  Is  used  in  the  BEM  result  shown  In 
Figure  6.2  so  that  results  from  both  methods  can  be  more  reasonably  com¬ 
pared.  Any  differences  resulting  from  the  use  of  either  a  single  point 
or  a  two-point  excitation  should  be  more  evident  In  the  results  with  tar¬ 
get  points  closer  to  the  source.  Figures  6.3  and  6.4  show  the  BEM 
results  for  such  close  range  targets  (d/X  from  0.5  to  1.0)  using  a 
two-point  and  a  single  point  excitation  respectively.  Only  a  close  look 
at  the  case  of  lowest  d/X  reveals  that  the  displacement  amplitude  In  the 
case  of  a  two-point  excitation  Is  slightly  lower.  The  difference  Is, 
however,  very  small  and  It  can  be  concluded  that,  for  d/X  of  0.5  and  lar- 


Figure  6.4  -  2-0  BEM  Free  Field,  Single  Point  Source,  D  *  0.0 


ger,  a  single  point  and  a  two-point  excitation  give  the  same  results. 
Therefore,  a  single  point  excitation  as  originally  formulated  will  be 
used  in  all  8EM  analyses  presented  later. 

Figure  6.5  shows  the  FEM  results  for  the  same  close  range  targets 
which  are  still  smoother  than  the  BEM  results  shown  In  Figures  6.3  and 
6.4.  The  Interesting  point  Is  that  the  displacement  amplitudes  in  this 
range  are  higher  than  the  corresponding  BEM  results. 

Figure  6.6  also  shows  the  close  range  results  from  the  BEM  using 
a  single  point  excitation  as  In  Figure  6.4,  but  with  soil  damping  ratio 
of  0.05.  The  peak  amplitudes  are  noticeably  lower  than  for  the  case 
without  damping  (Figure  6.4),  and  the  curves  are  also  smoother  similarly 
to  the  FEM  result. 

Figures  6.7,  6.8,  and  6.9  show  the  BEM  results  at  long  distance 
targets  for  the  case  of  soil  damping  ratios  D  of  0.0,  0.02,  and  0.05 
respectively.  By  comparing  Figure  6.7  to  Figure  6.2  at  the  same  d/X,  the 
Insignificant  effect  of  a  two-point  versus  a  single  point  excitation  is 
confirmed.  The  soil  damping  also  has  the  same  effects  of  the  results  of 
the  long  distance  targets  as  In  the  short  range,  I.e.,  the  amplitude  is 
lower  with  the  smoother  and  rounder  curves  resembling  the  FEM  results. 
This  makes  the  estimation  of  the  wave  arrival  time  more  difficult,  unlike 
the  case  without  soil  damping  In  Figure  6.7  where  the  time  of  arrival  of 
the  wave  Is  very  well  defined  by  the  sharp  rise  toward  the  first  positive 
peak.  From  all  results  above,  It  appears  that  the  time  of  arrival  of  the 
wave  decreases  with  Increasing  soil  damping. 


Using  log-log  scale,  the  amplitudes  of  the  first  positive  peak 
are  plotted  versus  d/X  for  different  soil  damping  ratios,  as  shown  In 
Figure  6.10.  The  amplitudes  from  the  FEM  results  (without  damping)  are 
also  plotted  together  In  the  same  figure.  For  the  BEM  results,  the  plot 
Is  a  straight  line  with  a  slope  equals  to  -1/2  when  there  Is  no  soil 
damping.  This  follows  the  geometrical  damping  law  for  the  cylindrical 
wavefront.  This  curve,  however,  deviates  from  a  straight  line  as  d/X 
decreases,  a  point  which  will  be  discussed  later.  In  the  cases  with  soil 
damping,  the  waves  attenuate  at  a  faster  rate  as  shown  In  the  plots  with 
the  lines  curving  downward.  As  the  distance  d  becomes  smaller,  there  are 
less  differences  In  the  peak  amplitudes  among  all  cases. 

Figure  6.10  also  reveals  another  distinct  character  of  a  cylin¬ 
drical  wavefront.  In  the  figure,  the  log-log  plots  of  the  attenuation 
for  all  soil  dampings  have  maximum  curvature  at  a  distance  of  approxi¬ 
mately  one  wavelength  (d/X  =  1).  As  the  distance  becomes  smaller,  all 
plots  show  a  reverse  curvature  and  tend  to  approach  each  other.  Further 
Investigations  reveal  that.  In  the  case  of  a  pure  shear  wave,  as  caused 
by  an  antiplane  excitation  (SH-wave),  there  Is  no  break  In  the  curvature 
and  the  plot  for  the  case  without  damping  Is  a  straight  line  throughout. 
The  deviation  from  a  straight  line  as  d/X  becomes  smaller  In  Figure  6.10 
Is,  therefore,  due  to  the  combination  of  SV-and  P-waves  caused  by  an 
Inplane  excitation. 

Figure  6.10  shows  further  that  the  wave  amplitudes  from  the  FEN 
solution  attenuate  faster  than  those  from  the  BB1  (both  have  no  damping). 
The  amplitude  from  the  FEM  Is  slightly  larger  for  d/X  below  2,  and 


becomes  smaller  than  the  amplitude  from  the  BEM  for  larger  distances. 
The  Important  point  to  note  Is  that  the  plot  of  the  FEM  does  not  show  any 
straight  line  portion  as  In  the  BEM  even  though  there  Is  no  soil  damping. 
From  these  results  together  with  the  curves  of  the  displacement  versus 
time  presented  earlier,  It  appears  that  the  finite  element  model  behaves 
like  a  system  with  a  small  amount  of  material  damping  on  the  order  of  1 
percent. 

6.1.2  ThrM-Dimmnsional  Case 

The  three-dimensional  free  field  displacements  for  different 
values  of  d/X  are  shown  In  Figures  6.11,  6.12,  and  6.13  for  D  -  0.0, 
0.02,  and  0.05  respectively.  All  results  show  the  early  motion  starting 
at  the  arrival  time  of  the  plane  P-wave  which  attenuates  very  fast  as  the 
distance  Increases.  In  the  case  D  =  0.0,  the  results  show  again  a 
well-defined  time  of  arrival  of  the  wave  while  In  the  case  with  soil 
damping,  the  results  become  smoother,  and  It  also  appears  that  the 
arrival  time  of  the  S-wave  decreases  with  increasing  soil  damping. 

While  the  points  mentioned  above  are  all  similar  to  those  of  the 
two-dimensional  case  shown  In  Figure  6.7  to  Figure  6.9,  the  main  differ¬ 
ences  between  the  two  cases  are  obvious.  In  the  three-dimensional  case, 
the  positive  and  the  negative  peaks  are  comparable  In  magnitude,  and  the 
displacement  after  a  complete  cycle  Is  perfectly  zero  Instead  of 
decreasing  asymptotically.  The  shape  of  the  displacement  curve  In  this 
case,  therefore,  retains  closely  the  sinusoidal  shape  of  the  excitation 
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except  for  the  early  motion  part  which  almost  disappears  at  larger  dis¬ 
tances. 

Figures  6.14  and  6.15  show  the  three-dimensional  free  field  dis¬ 
placements  at  close  range  targets  (d/X  from  0.5  to  1.0)  for  D  =  0.0  and 
0.05  respectively.  They  still  reproduce  a  better  sinusoidal  shape  of  the 
excitation  with  more  balanced  negative  and  positive  peaks  than  the 
two-dimensional  case  shown  In  Figures  6.4  and  6.5.  In  addition,  the 
change  In  the  direction  of  the  motion  seems  to  be  more  abrupt  In  the 
three-dimensional  case  as  shown  In  the  curve  as  a  sharp  rise  at  the 
S-wave  arrival  time  In  the  case  D  =  0.0.  All  results  from  both  two  and 
three-dimensional  cases  also  show  that  the  period  T  of  the  motion  slight¬ 
ly  Increases  with  the  damping  ratio  of  soil. 

Figure  6.16  shows  the  attenuation  curves  for  all  soil  damping  In 
log-log  scale.  Similarly  to  the  two-dimensional  case  In  Figure  6.10,  all 
results  have  the  maximum  curvature  at  d/X  =  1.  The  plots  show  a  reverse 
curvature  and  tend  to  approach  each  other  at  small  values  of  d/X  The 
only  difference  from  the  two-dimensional  case  Is  the  steeper  slope  which 
Indicates  the  faster  rate  of  attenuation  In  the  three-dimensional  case. 
For  D  =  0.0  and  d/X  larger  than  2,  the  plot  shows  a  straight  line  portion 
with  the  slope  equal  to  -1.  This  agrees  with  the  geometrical  damping  law 
for  a  spherical  wavefront. 


6.2  EFFECT  OF  AN  INCLUSION  ON  ATTENUATION 


Attenuation  studies  were  conducted  next  with  one  square  incla 
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sion  (cubic  In  the  three-dimensional  case)  having  side  dimension  'a1 
located  between  the  excitation  point  and  the  target.  Preliminary  ana¬ 
lyses  for  small  Inclusions  (a/X  up  to  1/4)  Indicated  that  the  results 
were  insensitive  to  the  position  of  the  Inclusion  between  the  excitation 
and  the  target  points.  No  studies  were  performed  for  the  case  of  larger 
Inclusions  where  the  position  of  the  Inclusion  may  have  some  effects  on 
the  results.  For  all  cases  to  be  consistent,  the  Inclusion  was  located 
at  a  fixed  position  varying  the  position  of  the  target  In  order  to  obtain 
results  for  different  values  of  d/X  from  a  single  run.  Various  sizes  of 
the  Inclusion  relative  to  the  wavelength  X  were  used.  The  Inclusion  had 
no  damping  and  Its  properties  (shown  with  subscript  'b')  relative  to  soil 
were 

csb/c$  =  10.0  ,  Pjj/p  *  1.0 
where  c$  Is  the  shear  wave  velocity,  and  p  Is  the  unit  mass. 


6.2.1  Two-Dimensional  Case 

The  displacements  for  the  case  of  one  square  Inclusion  with  size 
a/X  =  1/8,  and  soil  damping  ratios  D  of  0.0,  0.02,  and  0.05  are  shown 
In  Figures  6.17,  6.18,  and  6.19  respectively.  Only  a  close  comparison 
with  the  corresponding  free  field  case  can  reveal  a  slightly  smaller 
amplitude  In  this  case.  An  Inclusion  of  this  size  Is  too  small  In 
relation  to  the  wavelength  to  have  any  significant  effects.  This  can 
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also  be  seen  In  the  results  with  soil  damping  and  the  attenuation  curves 
shown  In  Figure  6.20  which  are  almost  the  same  as  In  the  free  field  case. 

When  the  size  of  the  Inclusion  Is  doubled  (a/X  =  1/4),  the  shape 
of  the  displacement  curves  Is  still  similar  to  those  of  the  free  field 
case  but  the  amplitudes  decrease.  This  can  be  seen  In  Figures  6.21, 
6.22,  and  6.23  for  damping  ratios  of  0.0,  0.02,  and  0.05  respectively. 
The  soil  damping  still  has  the  same  effects  of  smoothenlng  the  motion  and 
reducing  the  wave  arrival  time.  The  attenuation  curves  shown  in  Figure 
6.24  still  exhibit  the  same  characteristics  of  wave  attenuation  as  In  the 
free  field  except  that  the  amplitude  Is  lower. 

As  the  Inclusion  size  Is  further  Increased  (a/X  =  1/2),  the  over¬ 
all  amplitudes  are  reduced  even  more  as  shown  In  Figures  6.25,  6.26,  and 
6.27  for  D  =  0.0,  0.02,  and  0.05  respectively.  Figure  6.28  shows  that 
the  attenuation  curves  In  this  case  still  have  shapes  similar  to  those  of 
the  free  field  but  with  lower  amplitude. 

When  comparing  the  results  from  all  sizes  of  Inclusion,  It  can  be 
noticed  that  displacement  versus  time  curves  are  shifted  more  to  the  left 
as  the  Inclusion  getting  larger.  This  Indicates  a  faster  wave  velocity 
which  Is  reasonable  because,  with  a  larger  Inclusion,  the  wave  will  trav¬ 
el  less  In  soil  and  more  In  an  Inclusion  which  has  higher  shear  wave 
velocity.  The  effects  of  the  Inclusion  on  the  wave  arrival  time  will  be 
discussed  In  the  next  chapter. 

A  close  comparison  between  the  displacement  versus  time  curves 
from  each  Inclusion  size  also  reveals  that  the  period  T  of  the  motion 
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Figure  6.25  -  2-D  Displacement,  One  Inclusion,  a/x  =  1/2,  D  *  0.0 
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slightly  Increases  with  the  Inclusion  size  and  for  each  size,  this  period 
also  Increases  with  soil  damping  as  In  the  free  field. 

The  ratios  of  the  peak  amplitude  between  the  case  with  an  Inclu¬ 
sion  and  the  free  field  case  are  plotted  against  d/X  for  all  three  sizes 
of  Inclusion  (a/X  of  1/8,  1/4,  and  1/2).  These  plots  are  shown  In  Fig¬ 
ures  6.29,  6.30,  6.31,  and  6.32  for  the  damping  ratios  D  of  0.0,  0.02, 
0.05,  and  0.10  respectively.  All  plots  show  that  the  presence  of  an 
Inclusion  on  the  travelling  path  of  the  wave  will  reduce  the  wave  ampli¬ 
tude.  This  reduction  effect  Is  more  pronounced  with  larger  Inclusion 
size.  When  there  Is  no  soil  damping  the  reduction  In  amplitude  Is  almost 
constant  In  the  range  d/X  from  2  to  20.  These  reductions  are  about  35, 
15,  and  5  percent  for  a/X  of  1/2,  1/4,  and  1/8  respectively  as  shown  In 
Figure  6.29.  As  soil  damping  Increases,  the  amplitude  reduction 
decreases  especially  at  large  distance  as  shown  In  Figures  6.30  to  6.32. 
The  amplitude  ratios  A/A^  approach  the  same  value  as  d/X  Increases.  In 
all  cases,  the  maximum  reduction  In  amplitude  occurs  at  the  distance 
about  1.25  of  the  wavelength  X 


6.2.2  Three-Dimensional  Cat* 

The  displacement  versus  time  curves  for  the  case  with  one  cubic 
Inclusion  having  a  side  dimension  to  wavelength  ratio  a/X  =  1/8  are  shown 
in  Figures  6.33,  6.34,  and  6.35  for  damping  ratios  of  0.0,  0.02,  and  0.05 
respectively.  In  order  to  keep  the  computation  cost  within  a  reasonable 
limit,  the  truncation  of  the  transfer  function  was  made  at  a  slightly 
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Figure  6.31  -  Amplitude  Ratios  between  Cases  with  Inclusion  and 
Free  Field,  D  =  0.05,  2-D 
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lower  frequency  than  required  to  obtain  a  smooth  curve.  This  only  causes 
a  slight  Irregularity  on  some  portions  of  the  displacement  curve  In  the 
case  without  damping,  as  shown  In  Figure  6.33,  while  no  significant  accu¬ 
racy  Is  sacrificed.  In  the  case  with  damping,  most  of  the  higher  fre¬ 
quency  components  are,  however,  going  to  be  damped  out,  and  the  trun¬ 
cation  of  the  transfer  function  at  a  slightly  lower  frequency  has, 
therefore,  no  significant  effects  of  the  smoothness  of  the  results. 

Except  for  the  less  smooth  curve  due  to  the  truncation  effect  as 
mentioned,  there  are  no  significant  differences  between  the  results 
above  and  the  free  field  results  shown  In  Figures  6.11  to  6.13,  which 
Indicates  that  this  size  of  Inclusion  (a/X  =  1/8)  Is  too  small  to  have 
any  significant  effects.  The  attenuation  curves  for  this  case  (Figure 
6.36)  are,  therefore,  almost  Identical  to  those  of  the  free  field  (Figure 
6.16). 

When  the  Inclusion  size  Is  doubled  (a/X  =  1/4),  the  peak  ampli¬ 
tude  decreases  slightly  as  shown  In  Figures  6.37  and  6.38  for  the  case 
with  damping  ratios  of  0.0  and  0.02  respectively.  With  damping  ratios  of 
0.05,  shown  In  Figure  6.39,  the  amplitude  decrease  can  be  detected  only 
at  close  distances  (d/X  of  4  to  8)  beyond  which  the  amplitude  are  all  the 
same  with  the  case  a/X  =  1/8  and  the  free  field  case.  This  Indicates  a 
larger  effect  of  the  Inclusion  with  closer  targets  and  less  soil  damping. 
Figure  6.40  shows  the  corresponding  attenuation  curve  which  Is  still 
basically  a  straight  line  In  the  case  D  =  0.0. 

In  the  results  shown  In  Figure  6.41,  the  size  of  the  Inclusion  Is 
further  Increased  (a/X  =  1/2).  Because  of  the  limit  In  the  computer  sto- 
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igure  6.39  -  3-D  Displacement,  One  Inclusion,  a/A  =1/4,  D  =  0.05 
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rage  available,  the  number  of  degrees  of  freedom  had  of  be  kept  small  by 
using  larger  elements  which  are  below  the  limit  required  for  accuracy. 
The  side  dimension  of  the  boundary  element  used  in  this  case  Is  X/4 
Instead  of  X/8  as  usual.  The  results,  as  expected,  are  not  as  accurate 
especially  when  the  target  point  Is  close  to  the  inclusion  as  shown  In 
the  curve  when  d/X  =  4.  Despite  the  loss  In  accuracy,  the  results  still 
Indicate  a  larger  amplitude  reduction  in  the  closer  range  target. 

The  ratios  of  the  peak  amplitude  between  the  case  with  an  inclu¬ 
sion  and  the  free  field  case  are  plotted  against  d/X  In  the  same  manner 
as  In  the  two-dimensional  case  for  a/X  of  1/8  and  1/4.  These  results  are 
shown  In  Figures  6.42  to  6.45  for  damping  ratios  of  0.0  to  0.10  respec¬ 
tively.  For  a/X  =  1/8,  there  Is  no  reduction  in  amplitude  except  in  the 
case  without  damping  which  shows  a  decrease  In  amplitude  of  about  2  per¬ 
cent  at  d/X  =  4,  down  to  1  percent  at  d/X  =20.  With  larger  inclusion  size 
(a/X  =  1/4),  the  reduction  Is  more  evident.  This  reduction,  however, 
decreases  with  the  larger  distance  d/X  and  larger  damping  ratio  D.  With 
the  damping  ratio  D  =  0.10,  there  Is  no  reduction  In  amplitude  at  all 
distance  ranges  studied  as  shown  In  Figure  6.45. 

The  above  results,  as  compared  to  the  two-dimensional  case, 
reveals  that  the  amplitude  reduction  Is  caused  more  by  the  soil  damping 
than  by  the  presence  of  an  Inclusion.  This  explains  why  there  are  no 
differences  In  the  results  between  the  free  field  and  the  case  with 
Inclusion  when  the  damping  ratio  Increases. 


Figure  6.43  -  Amplitude  Ratios  between  Cases  with  Inclusion  and 
Free  Field.  0  *  0.02.  3-D 
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igure  6.44  -  Amplitude  Ratios  between  Cases  with  Inclusion  and 
Free  Field,  0  *  0.05,  3-D 
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Figure  6.45  -  Amplitude  Ratios  between  Cases  with  Inclusion  and 
Free  Field,  D  *  0.10,  3-D 
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6.3  EFFECT  OF  THE  MASS  AND  STIFFNESS  OF  AN  INCLUSION 

Only  the  stiffness  and  Inertia  of  the  Inclusion  which  are  the 
essential  system  properties  for  the  transmission  of  wave  motion,  are 
Investigated  In  this  section.  The  Inclusion  considered  here  is  a  square 
with  the  size  parameter  a/X  =  1/8.  Two  sets  of  results  are  presented. 
In  the  first  set,  the  Inclusion  Is  100  times  stlffer  than  the  medium 
(cs^/cs  =  10),  while  In  the  second  set,  the  Inclusion  Is  softer  than  the 
medium  (csl)/cs  =  0.5).  In  each  set,  the  unit  mass  ratios  between  the 
Inclusion  and  the  medium  are  0.5,  1.0,  and  5.0. 

Figure  6.46  shows  the  plots  of  amplitude  versus  distance  In 
log-log  scale  for  the  case  of  the  stiff  Inclusion.  The  plots  show  a 
slightly  higher  amplitude  with  the  higher  unit  mass  ratio.  The  rate  of 
amplitude  attenuation  Is,  however,  the  same. 

In  the  case  of  an  Inclusion  that  Is  softer  than  the  medium,  the 
effect  of  the  Inclusion  unit  mass  on  the  amplitude  and  the  attenuation 
characteristics  are  the  same  as  In  the  case  of  a  stiff  Inclusion;  but, 
the  Increase  In  amplitude  with  the  unit  mass  ratio  Is  slightly  larger  as 
can  be  seen  In  Figure  6.47.  By  comparing  this  figure  to  Figure  6.46,  It 
can  be  noticed  that,  for  the  same  unit  mass  ratio,  the  case  of  a  stlffer 
Inclusion  shows  a  slightly  smaller  amplitude  than  the  soft  Inclusion. 
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CHAPTER  7 

EFFECT  OF  INCLUSIONS  ON  ARRIVAL  TIME 

Stress  wave  propagation!  techniques  are  often  used  In  geotechni¬ 
cal  engineering  for  the  exploration  of  soil  deposits  and  determining 
soil  moduli  In  the  field  from  measurements  of  wave  velocities.  In  a  non- 
homogeneous  medium  or  a  medium  with  inclusions  having  different  proper¬ 
ties,  such  as  an  alluvial  deposit,  the  velocity  of  the  propagating  wave 
would  be  different  depending  on  a  complicated  dynamic  Interaction 
between  the  medium  and  the  Inclusions.  This  Interaction  will  be  a  func¬ 
tion  of  the  physical  properties,  shape,  and  size  of  the  Inclusions  as 
well  as  their  number  and  their  arrangement  In  the  cluster.  As  mentioned 
at  the  beginning  of  Chapter  6,  extensive  studies  on  these  factors  Involve 
a  considerable  expense  and  time.  Consequently,  only  some  major  cases  for 
the  two-dimensional  model  will  be  Investigated  and  presented  In  this 
chapter. 


7.1  INTERPRETATION  AND  PRESENTATION  OF  THE  ARRIVAL  TIME 

Since  the  arrival  time  Is  an  Interpreted  value  from  the  plots  of 
motion  (which  may  be  the  displacement,  velocity,  or  acceleration)  versus 
time,  the  Interpretation  method  has  to  be  consistent  In  order  to  have 


reliable  results.  Although  the  determination  of  the  wave  propagation 
velocity  from  the  Interval  travel  time  between  two  target  points  proves 
to  be  somewhat  easier  and  more  reliable  than  from  the  first  arrival  time 
at  one  location,  the  latter  approach  will  be  used  here  since  the  main 
point  of  interest  In  this  work  Is  the  relative  change  of  the  wave  veloci¬ 
ty  (or  arrival  time)  with  the  presence  of  Inclusions.  As  long  as  the 
method  Is  consistent,  the  results  should  show  the  same  Inclusion  effect. 

The  determination  of  the  first  arrival  time  of  the  shear  wave 
using  the  displacement  versus  time  trace  for  the  case  of  a  finite  domain 
with  a  point  force  excitation  (Suddhlprakarn,  1983)  showed  that  the 
first  arrival  time  could  be  either  selected  as  the  zero  crossing  point 
(when  the  positive  motion  starts),  or  the  point  where  the  major  upward 
movement  starts  (as  often  used  In  experimental  work).  The  differences  in 
the  shear  wave  velocity  between  the  above  two  criteria  Is  about  five  per¬ 
cent.  To  justify  the  accuracy  of  these  two  approaches  for  the  case  of  an 
Infinite  medium  using  the  boundary  element  method,  the  free  field  dis¬ 
placement  for  short  range  and  long  range  targets  are  plotted  together  In 
Figure  7.1  and  Figure  7.2  respectively.  The  dimensionless  displacement 
and  time  defined  In  Chapter  6  are  used,  which  leads  to  a  different  hori¬ 
zontal  scale  for  curves  of  different  targets  depending  on  the  tar¬ 
get-source  distance.  The  theoretical  arrival  time  of  the  shear  wave  to 
all  targets  In  the  free  field,  however,  will  always  correspond  to  a  value 
of  c$t/d  of  one. 

The  results  for  the  close  range  targets  In  Figure  7.1,  whose  hor¬ 
izontal  scale  Is  highly  enlarged  show  that  the  zero  crossing  point  does 


not  represent  the  first  arrival  time  as  well  as  the  point  where  the  major 
upward  motion  starts.  This  point  can  be  located  at  the  lowest  point  on 
the  early  motion  part.  As  the  distance  from  the  source  Increases,  the 
early  motion  starts  to  show  a  double  curvature  trace  with  a  minimum  point 
ahead  of  the  theoretical  arrival  time  as  shown  In  Figure  7.2.  In  this 
case,  It  Is  clear  that  the  arrival  time  would  be  best  selected  at  the 
point  where  the  big  upward  movement  starts. 

Both  results  for  the  close  range  and  the  long  range  targets, 
therefore,  reveal  that  the  point  where  the  major  upward  motion  starts  Is 
always  the  best  point  to  represent  the  first  shear  wave  arrival  time 
Irrespective  of  the  shape  of  the  early  motion  curve.  This  approach  will, 
therefore,  be  used  throughout  this  work. 


7.2  EFFECT  OF  THE  INCLUSION  SIZE  AND  WAVELENGTH 

Some  results  in  Chapter  6  had  already  shown  a  trend  of  a  faster 
shear  wave  velocity  with  the  Increase  In  the  size  of  the  Inclusion.  This 
Is  reasonable  because  the  wave  travel  more  In  the  Inclusion  which  has  a 
higher  wave  velocity.  On  the  other  hand.  If  the  Inclusion  has  a  lower 
shear  wave  velocity,  the  effect  should  reverse.  It  Is  obvious  that  if 
the  travelling  path  of  the  wave  In  an  Inclusion  "d^"  makes  up  a  large 
portion  of  the  total  travelling  path  "d",  the  effect  of  the  inclusion 
size  would  be  much  more  pronounced  as  shown  In  Figure  7.3.  In  this  fig¬ 
ure,  the  displacements  for  square  Inclusions  of  different  sizes  (a/X  = 
1/4  and  1/8)  are  plotted  together  with  the  free  field  displacement,  all 


One  Inclusion 


csb/cs 

d/x 


1 


=  10 
=  0.5 


i 


at  d/X  =  0.5.  The  shear  wave  velocity  of  the  Inclusions  In  this  case  Is 
ten  times  faster  than  that  of  the  medium  while  other  properties  are  the 
same.  The  figure  clearly  shows  that  the  displacement  curve  shifts  to  the 
left  as  the  Inclusion  size  and  the  travelling  path  d^  Increase. 

Figure  7.4  shows  the  variation  of  the  dimensionless  variable 
c$t/d  versus  d^/d  (d^  length  of  the  Inclusion  In  the  direction  of  propa¬ 
gation,  d  distance  between  the  source  and  the  target)  for  a 
one-dimensional  situation  where  only  one  ray  would  travel  through  the 
distance  d-db  with  a  velocity  c$  and  through  the  distance  dfa  with  a 
velocity  c$b.  For  this  simplified  situation,  one  can  define  an  equiv¬ 
alent  distance 
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When  the  Inclusion  Is  stiffen  than  the  surrounding  medium  and 
therefore  csb/c$  Is  larger  than  one,  this  formula  should  provide  a  lower 
bound  to  the  solution  for  a  two  dimensional  situation.  On  the  other 
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hand,  when  the  Inclusion  Is  softer  than  the  medium  the  formula  should 
provide  an  upper  bound.  Curves  based  on  this  formula  will  be  used  for 
the  purpose  of  comparison  when  studying  the  effect  of  the  Inclusion  on 
the  time  of  arrival  or  apparent  wave  velocity  for  different  Inclusions. 

The  effect  of  the  inclusion  should^ b*  larger  when  dealing  with 
ratios  of  d^/d  close  td  unity.  This  implies  that  both  a  source  and  a 
target  point  are  relatively  close  to  the  inclusion.  In  this  case,  more 
boundary  elements  (or  smaller  element  sizes)  may  be  necessary  to  main¬ 
tain  the  same  level  of  accuracy.  Plots  of  the  arrival  time  versus  the 
ratio  d^/d  for  the  case  a/X  =  1/16  using  three  different  element  sizes 
are  shown  in  Figure  7.5.  The  case  of  the  element  length  l  =  X/16  corre¬ 
sponds  to  a  single  element  per  face  of  the  Inclusion,  that  Is  to  say  a 
total  of  four  boundary  elements.  The  case  l  =  X/32  corresponds  to  two 
elements  per  face  and  l  =  X/64  to  four  elements  per  face  (  a  total  of  16 
elements).  It  can  be  seen  that,  as  the  element  size  decreases,  the 
arrival  time  curve  shifts  up  and  down  (there  Is  no  monotonic  conver¬ 
gence),  and  to  reach  the  final  position,  the  element  size  has  to  be  very 
small.  In  our  case,  this  final  position  was  determined  by  judging  from 
the  results  using  two  or  three  different  element  sizes. 

The  first  arrival  times  cst/d  for  different  values  of  the  ratio 
d/d^  are  shown  In  Table  7.1  for  three  different  wacelengths  (a/X  =  1/16, 
1/8,  and  1/4).  The  values  from  this  table  are  plotted  In  Figure  7.6 
together  with  the  curve  for  the  1-D  ray  for  comparison.  These  plots  show 
that  the  arrival  time  Is  larger  as  expected  for  the  longer  wavelength. 
This  Is  because  as  the  wavelength  Increases  In  relative  to  the  inclusion 
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Values  of  cst/d  at  the  First  Arrival  of  the  Shear  Wave 


N.  a/x 

d/db\ 

1/4 

1/8 

1/16 

1.75 

0.54 

0.56 

0.57 

2.00 

0.62 

0.64 

0.65 

2.50 

0.75 

3.00 

0.79 

0.82 

0.85 

4.00 

0.90 

5.00 

0.88 

0.93 

0.94 

6.00 

0.94 

7.00 

0.90 

0.94 

0.96 

db=a  (square  inclusion) 


size,  most  of  the  wave  will  propagate  through  the  medium  Instead  of  the 
Inclusion  and  the  effect  of  the  Inclusion  Is  thus  smaller.  As  a  result, 
the  arrival  time  curve  tends  to  be  nearer  to  the  upper  bound  at  c$t/d  of 
one  (for  the  free  field  case).  The  results  show  that  the  decrease  of  the 
wavelength  a/X  from  1/16  to  1/8  does  not  change  the  result  as  much  as  the 
decrease  of  a/X  from  1/8  to  1/4.  This  Indicates  the  smaller  effect  from 
the  inclusion  when  a/X  Is  less  than  1/8. 

7.3  EFFECT  OF  INCLUSION  SHAPE 

In  the  study  of  the  effect  of  the  shape  of  the  Inclusion,  It  Is 
more  difficult  to  standardize  the  procedure  because  of  many  types  of 
shapes  that  can  be  considered.  In  addition,  different  Inclusion  shapes 
mean  also  different  sizes.  In  an  attempt  to  neutralize  the  effect  of 
wavelength  and  the  Inclusion  size  as  much  as  possible,  two  groups  of 
shapes  considered  In  this  study  are  selected  on  the  basis  of  two  differ¬ 
ent  criteria: 

1.  All  Inclusions  have  the  same  length  of  the  wave  path  measured  along 
the  direction  of  propagation,  (which  Is  a  one-dimensional  consider¬ 
ation).  The  Inclusion  shapes  In  this  group  are  a  square,  a  vertical 
rectangle,  a  horizontal  rectangle  (having  the  sane  proportion  as  the 
vertical  one)  which  Is  considerably  smaller  In  order  to  meet  the  cri¬ 
terion,  and  lastly,  a  long  vertical  rectangle  with  a  height  of  1.125 
of  the  wavelength. 


150 


2.  All  Inclusions  have  the  same  cross-sectional  area  (which  Is  more  of  a 
two-dimensional  consideration).  The  Inclusion  shapes  In  this  group 
are  a  vertical  rectangle,  a  horizontal  rectangle,  a  rhombus,  a 
square,  and  a  circle  which  Is  actually  a  dodecagon  (12  sides  polygon) 
because  It  Is  approximated  by  12  straight  elements. 

The  shapes  and  sizes  of  the  inclusion  of  these  two  sets  with  the 
relative  size  of  the  wavelength  X  are  shown  In  Figure  7.7.  As  in  the 
last  section,  the  shear  wave  velocity  of  the  Inclusion  used  In  this  study 
Is  ten  times  that  of  the  medium,  while  another  properties  are  the  same, 
and  there  Is  no  damping  In  any  of  the  materials. 

The  arrival  time  versus  d/db  curves  for  the  Inclusions  under  the 
first  set  are  shown  in  Figure  7.8.  The  arrival  time  Is  largest  (smallest 
apparent  wave  velocity,  closest  to  that  of  the  medium)  for  the  case  of  a 
horizontal  rectangle  whose  size  Is  much  smaller  than  the  rest  In  order  to 
satisfy  the  criterion.  The  wave  velocity  Increases  In  the  case  of  a 
square,  a  vertical  rectangular,  and  a  long  vertical  rectangular  Inclu¬ 
sions,  consecutively.  At  larger  distance,  there  Is  less  difference  In 
the  wave  velocity  between  the  case  of  a  normal  length  and  a  long  vertical 
rectangle.  These  results  Indicate  that  the  dimension  of  the  Inclusion 
perpendicular  to  the  direction  of  wave  propagation  affects  also  the  sol¬ 
ution  for  the  two-dimensional  case.  As  this  dimension  Increases,  the 
solution  should  approach  the  one-dlmenslonal  curve. 

Figure  7.9  shows  the  arrival  time  curves  for  the  Inclusions  In 
the  second  set.  It  can  be  seen  that  the  wave  velocity  Is  slowest  for  a 
rhombus  and  fastest  for  a  square  while  the  case  of  a  circular  Inclusion 
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Is  In  between.  The  case  of  a  horizontal  rectangle  shows  about  the  same 
wave  velocity  as  that  of  the  circle.  The  curve  for  the  vertical  rectan¬ 
gle  shows  a  much  faster  wave  velocity  at  small  d/dfa  (close  to  the  case  of 
a  square  Inclusion)  and  slower  as  d/d^  Increases. 

It  should  be  noted  that  of  all  five  Inclusion  shapes  considered, 
the  rhombus  and  the  circle  have  about  the  same  length  of  the  principal 
wave  path  d^,  both  of  which  are  slightly  longer  than  that  of  a  square 
Inclusion.  These  three  Inclusions  are,  however,  comparable  In  overall 
proportions.  On  the  other  hand,  the  vertical  and  the  horizontal  rectan¬ 
gles,  each  have  significantly  different  path  length  from  the  rest. 
Because  of  this,  the  results  from  both  rectangular  shapes  cannot  compare 
well  to  the  other  cases  since  all  plots  are  normalized  on  the  basis  of 
the  value  of  the  wave  path's  length. 

The  value  of  c$t/d  at  the  first  arrival  time  of  the  shear  wave  at 
different  distances  from  source  (normalized  with  the  fixed  wavelength  X) 
are  tabulated  In  Table  7.2  and  Table  7.3. 


7.4  MULTIPLE  INCLUSION 

As  In  the  case  of  the  Inclusion  shape,  It  Is  difficult  to  select 
meaningful  parameters  to  use  as  a  basis  for  comparison  between  different 
arrangements  of  the  Inclusion  cluster.  The  basis  used  In  this  section, 
which  may  not  be  the  best.  Is  that  the  total  distance  In  the  Inclusions 
(db)  along  the  principal  axis  of  wave  propagation  Is  the  same,  as  In  the 
first  criteria  In  the  last  section.  In  this  section  all  Inclusions  con- 


Table  7.2  -  First  Arrival  Time  c  t/d  For  Inclusions  Under  Shape 
Criterion  1  s 
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sldered  are  either  rectangle  or  square  with  the  least  dimension  "a"  equal 
to  1/8  of  the  wavelength.  The  properties  of  the  medium  and  the  Inclu¬ 
sions  are  as  In  the  last  section. 

Figure  7.10  compares  the  arrival  time  curves  between  the  case  of 
a  vertical  rectangular  Inclusion  and  a  group  of  three  square  Inclusions 
aligned  perpendicularly  to  the  principal  wave  path.  In  both  cases,  the 
total  area  of  the  Inclusions  are  also  the  same.  The  previous  results  for 
the  case  of  one  square  Inclusion  are  also  shown  for  comparison.  Of  all 
three  cases,  the  wave  Is  fastest  In  the  case  of  a  rectangular  Inclusion 
and  slowest  In  the  case  of  one  square  Inclusion.  The  results  Indicate 
that  the  effect  on  the  wave  velocity  Is  maximum  when  the  inclusion  Is  in 
one  large  continuous  piece  serving  as  a  long  barrier  so  that  most  of  the 
wave  propagates  through  the  Inclusion.  It  appears  that,  for  the  case  of 
an  Inclusion  group,  some  of  the  wave  can  propagate  through  the  medium 
between  Inclusions  causing  a  smaller  effect  on  the  wave  velocity. 

In  Figure  7.11,  a  similar  comparison  Is  made  between  the  case  of 
the  horizontal  rectangular  Inclusion  and  a  group  of  three  square  Inclu¬ 
sions  oriented  along  the  principal  path  of  the  propagating  wave.  Both 
cases  have  the  same  total  length  of  wave  path  In  the  Inclusions.  The 
results  show  the  faster  wave  velocity  In  the  case  of  a  continuous  rectan¬ 
gular  Inclusion  as  In  the  previous  case. 

Figure  7.11  shows  also  the  arrival  time  curve  for  a  group  of  nine 
square  Inclusions  arranged  In  a  3x3  block,  as  shown.  By  comparing  the 
results  to  those  for  a  group  of  three  inclusions,  It  can  be  seen  that  the 
addition  of  the  extra  rows  of  Inclusions  Increases  the  wave  velocity 
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Figure  7.11  -  First  Arrival  Time  for  the  Inclusion  Cluster  vs.  Single  fnclusion 


slightly,  but  It  Is  still  slower  than  from  the  case  of  one  continuous 
Inclusion. 

As  In  the  case  of  the  vertical  and  the  horizontal  rectangle,  the 
results  between  the  inclusion  cluster  that  are  arranged  vertically  and 
horizontally  cannot  be  compared  very  well  as  shown  In  the  figures 
because  of  the  different  length  of  the  total  principal  wave  path  "d^"  in 
the  Inclusions  and,  consequently,  the  dimensionless  parameter  d/dfa.  On 
the  basis  of  the  same  source-target  distance  "d" ,  the  horizontally 
aligned  group  has  more  effect  of  the  wave  velocity  due  to  the  longer  "d^" 
This  can  be  seen  In  Table  7.4  and  Table  7.5  (at  d/X  =  0.75)  where  the  val¬ 
ues  of  c$t/d  at  the  first  arrival  time  of  the  shear  wave  at  different 
distances  (relative  to  the  fixed  wavelength  X)  are  tabulated  for  differ¬ 
ent  Inclusion  clusters.  The  values  for  the  single  inclusions  are  also 
shown  for  comparison. 

It  can,  therefore,  be  concluded  that  Increasing  the  area  of 
Inclusions,  either  as  a  larger  Inclusion  (compared  with  the  wavelength) 
or  a  larger  number  of  inclusions  In  a  group,  will  Increase  the  effect  on 
the  wave  velocity.  More  Importantly,  for  the  equivalent  total  Inclusion 
area,  one  large  continuous  Inclusion  will  affect  the  wave  velocity  more 
than  a  group  of  several  smaller  Inclusions  (with  the  same  total  area  and 
travelling  path). 


7.5  INCLUSION  STIFFNESS 


The  results  presented  earlier  dealt  only  with  inclusion  with  a 
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shear  wave  velocity  ten  times  (or  shear  stiffness  of  100  times)  that  of 
the  medium.  In  many  real  situations,  the  Inclusions  may  have  a  wide 
range  of  shear  wave  velocities  different  from  that  of  the  medium.  They 
may  represent  a  rock  or  a  boulder,  a  pocket  of  softer  material,  sand  or 
silt  In  a  soil  deposit  or  even  an  underground  structure.  This  section 
will  shows  the  results  of  the  shear  wave  arrival  time  for  the  case  of  a 
square  Inclusion  (a/X  =  1/8)  with  different  shear  wave  velocity  ratios 
c^/c  (0.1,  0.5,  5,  10,  and  20)  while  other  properties  of  the  Inclusion 
and  the  medium  are  the  same.  Both  materials  have  no  damping. 

In  the  previous  studies  with  an  Inclusion  having  csb/cs  of  10* 
the  element  size  used,  which  Is  1/8  of  the  wavelength,  Is  normally  ade¬ 
quate  for  an  accurate  result  except  when  the  source  and  the  target  are 
very  close  to  the  Inclusion.  In  the  study  of  Inclusions  with  different 
stiffness.  It  was  also  found  that  the  element  size  1/8  of  the  wavelength 
works  well  until  the  ratio  csjJ/cs  goes  lower  than  0.5  where  the  results 
start  to  be  less  reliable.  This  is  demonstrated  In  Figure  7.12  where  the 
displacement  curves  for  an  Inclusion  with  csb/cs  of  0.1  using  two  differ¬ 
ent  element  sizes  are  plotted  together.  This  Is  because  the  surface  of 
the  less  stiff  Inclusion  cannot  retain  the  shape  as  In  the  case  of  the 
stlffer  Inclusion.  The  displacement  at  the  surface  Is,  therefore,  larg¬ 
er  and  less  uniform  and  cannot  be  accurately  represented  by  just  a  few 
elements.  In  this  section,  more  elements  are,  therefore,  used. 

Figure  7.13  shows  the  arrival  time  curves  for  different  Inclu¬ 
sion  stiffness  (Cjjj/Cj  of  0.1,  0.5,  5,  10,  and  20).  The  curve  for  the 
free  field  case  (csi/cs  =  1-0)  represented  by  a  horizontal  line  at 
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Figure  7.12  -  Effect  of  Element  Size  on  the  Displacement  in  the 
Case  of  Soft  Inclusion 
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cst/d  of  one.  The  vertical  scale  In  this  figure  Is  also  displaced 
slightly  to  accommodate  the  plots  for  the  case  of  c^/c  less  than  one. 
The  results  show  that  the  wave  velocity  Is  faster  with  the  Increase  In 
the  Inclusion  stiffness  as  expected.  The  arrival  times  for  the  cases 
where  c  ^/cs  equals  to  5,  10,  and  20  are  not  much  different  from  each 
other  and  these  curves  look,  similar  to  the  curves  of  a  1-D  ray  shown  In 
Figure  7.4,  but  somewhat  higher.  When  the  Inclusion  Is  less  stiff  than 
the  medium,  the  wave  velocity  reduces  slightly,  not  as  much  as  In  the 
case  of  a  1-0  ray.  This  can  be  seen  by  comparing  the  arrival  time  curve 
for  the  case  c  ./c$  =  0.1  In  Figure  7.13  and  Figure  7.4.  This  result, 
therefore,  reveals  that  the  softer  Inclusions  may  not  have  a  strong 
effect  In  reducing  the  overall  wave  velocity  since  the  wave  can  arrive 
faster  following  a  route  of  higher  wave  velocity. 


7.6  INCLUSION  MASS 

Most  of  the  Inclusions  considered  so  far  have  the  same  unit  mass 
as  the  surrounding  medium.  Actually,  the  unit  mass  of  the  Inclusion  can 
also  have  a  wide  range  of  values  from  as  low  as  an  air  void  to  a  heavy 
burled  structures. 

To  study  the  effect  of  the  unit  mass  of  the  Inclusion,  stiff  and 
soft  Inclusions  are  considered  as  In  Section  6.3.  The  stiff  and  soft 
Inclusions  considered  have  c^/Cj  of  10  and  0.5  respectively.  In  each 
case,  the  unit  mass  ratios  between  the  Inclusion  and  the  medium  are  0.5, 
1.0,  and  5.0. 
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Figure  7.14  shows  the  displacement  versus  time  plots  from  the 
case  of  a  stiff  Inclusion  for  all  three  unit  mass  ratios  mentioned  above. 
These  results  were  selected  from  the  target  at  d/X  =  0.625  because  at 
such  distance.  It  Is  not  too  close  to  an  Inclusion  surface  to  make  the 
results  lack  In  accuracy,  and  It  Is  yet  small  enough  that  the  change  of 
the  arrival  time  can  still  be  noticed  easily.  The  plots  show  that  the 
difference  In  the  arrival  times  of  the  shear  wave  between  the  unit  mass 
ratio  of  0.5  and  1.0  are  Insignificant.  With  a  further  Increase  of  the 
unit  mass  ratio  up  to  5.0,  the  Increase  In  the  arrival  time  Is  now  more 
obvious.  Figure  7.15  shows  the  same  type  of  plots  for  the  case  of  a  soft 
Inclusion  with  the  same  trends  of  results. 


Figure  7.15  -  Effect  of  the  Unit  Mass  of  the  Inclusion,  Soft 
Inclusion 
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CHAPTER  8 


CONCLUSIONS 


This  research  resulted  as  an  ex'anslon  of  analyses  performed  to 
support  an  experimental  Investigation  on  the  effect  of  the  state  of 
stresses  on  wave  velocities  In  soil  (Stokoe  et  al ,  1980).  In  the  ori¬ 
ginal  Investigation,  a  model  based  on  the  finite  element  method  (FEM) 
with  an  explicit  Integration  scheme  was  Implemented  to  simulate  the  pro¬ 
pagation  of  waves  Inside  a  big  trlaxlal  cube  containing  soil,  and  the 
Instrument  packages  as  rigid  Inclusions.  The  procedure  was  suitable  for 
that  particular  problem  because  of  the  finite  domain  and  the  possibility 
to  study  c  fferent  boundary  effects. 

In  this  work,  the  more  general  area  of  wave  propagation  In  a 
medium  with  Inclusions  having  different  properties  Is  addressed  by  con¬ 
sidering  a  full  space  so  that  boundary  effects  can  be  eliminated.  In 
this  case,  an  analysis  In  the  frequency  domain  using  the  boundary  element 
method  has  obvious  advantages  while  the  finite  element  model  from  the 
previous  work  was  modified  to  approximate  an  Infinite  medium  In  order  to 
compare  the  results  to  those  of  the  BEM. 

The  FEM  used  In  this  work  Is  based  on  the  displacement  method 
with  a  linear  basis  function  and  four-node  rectangular  elements.  The 
equations  of  equilibrium  for  a  finite  element  system  In  motion  were 
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derived  by  applying  Hamiltlon's  principle.  With  a  lumped  mass  matrix, 
neglecting  the  velocity  dependent  damping  forces  and  using  an  explicit 
Integration  scheme,  the  solution  can  be  marched  out  In  the  time  domain, 
at  the  element  level,  without  assembling  the  global  stiffness  matrix. 
The  time  step  of  Integration,  however,  has  to  be  sufficiently  small  to 
guarantee  stability.  In  this  FEM  model,  an  Inclusion  Is  represented  by  a 
single  element,  and  an  Infinite  medium  Is  approximated  by  extending  the 
boundaries  far  enough  so  that  they  have  Insignificant  effects. 

In  the  BEM,  Mie  formulation  for  both  two  and  three  dimensional 
cases  was  carried  out  In  the  frequency  domain.  Starting  with  the  Bet¬ 
ti -Maxwell  reciprocity  relation  and  a  fundamental  solution  In  transform 
space  for  a  unit  harmonic  point  force  In  an  Infinite  medium,  the 
expression  for  displacements  (so  called  Somlgllana's  Identity)  Is 
derived.  The  Integral  equation  Is  then  formed  by  taking  the  limiting 
case  of  Somlgllana's  Identity  for  the  point  on  the  boundary.  This  proce¬ 
dure  Is  carried  out  both  for  the  Inclusion  and  for  the  medium  with  cavi¬ 
ties  and  the  excitation  force.  Both  parts  are  then  combined  to  satisfy 
equilibrium  and  compatibility  conditions  leading  to  a  system  of 
equations  that  can  be  solved  for  displacements  and  tractions  at  the  boun¬ 
dary.  The  displacements  at  any  points  In  the  domain  can  be  solved  using 
again  Somlgllana's  Identity.  These  displacements  or  transfer  functions 
are  solved  at  different  frequencies  and  multiplied  by  the  Fourier  trans¬ 
form  of  the  loading  function.  The  Inverse  Fourier  transforms  of  this 
product  are  then  obtained  to  get  the  time  responses. 
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A  discrete  Fourier  transform  based  on  the  Fast  Fourier  transform 
(FFT)  algorithm  of  Cooley  and  Tukey  (1965),  is  used  In  this  work..  Since 
the  transfer  function  has  to  be  supplied  at  every  frequency  component  in 
the  FFT  and  the  exact  calculation  for  all  of  them  will  be  prohibitively 
expensive,  part  of  the  transfer  function  will  either  be  defined  by 
Interpolation  or  be  neglected  (truncation).  Preliminary  studies 
revealed  that  the  use  of  extrapolated  values  worsened  the  results.  Trun¬ 
cation  was,  therefore,  used  but  the  frequency  beyond  which  the  transfer 
function  Is  neglected  must  be  sufficiently  high  to  reproduce  the  excita¬ 
tion  frequency.  In  the  preliminary  studies,  other  parameters  involved 
In  frequency  domain  analyses,  and  the  size  of  the  finite  element  domain 
necessary  to  represent  an  Infinite  medium  were  also  investigated. 

The  main  body  of  this  study  focuses  on  the  effect  of  Inclusions 
on  the  wave  attenuation,  wave  form,  and  wave  velocity  (In  terms  of 
arrival  time)  for  shear  waves.  The  results  from  the  Investigation  lead 
to  the  following  conclusions: 

2-D  Free  Field  Attenuation:  BEM  versus  FEM 

1.  With  (a  shear  type  excitation  over  a  small  area,  which  can  be  consid¬ 
ered!  as  a  point  (as  used  throughout  this  work),  the  displacement 
curves  of  both  FEM  and  BEM  show  a  small  motion  starting  at  the 
arrival  time  of  the  P-wave.  This  early  motion  flattens  out  very  fast 
as  the  distance  from  the  source  Increases.  The  displacement  curve 
from  the  FEM  Is  smoother  and  decays  with  an  oscillatory  motion  after 


173 


the  loading  cycle  in  contrast  with  the  asymptotic  decay  obtained 
with  the  BEM.  Both  positive  and  negative  peaks  in  the  FEM  results 
are  almost  the  same  in  amplitude  while  for  the  BEM,  the  positive 
amplitude  is  considerably  larger  than  the  negative  one. 

2.  For  a  material  without  Internal  damping,  the  positive  peak  amplitude 
from  the  FEM  Is  larger  at  close  range  and  smaller  as  d/X  Increases 
than  that  obtained  from  the  BEM,  which  implies  that  the  attenuation 
rate  is  slightly  faster  for  the  former.  Furthermore,  the  plot  of 
amplitude  versus  distance  for  the  FEM  Is  a  curve  throughout.  This 
behaviour  looks  similar  to  that  of  systems  with  a  small  amount  of 
material  damping,  indicating  that  the  FEM  model  used  Introduces  some 
fictitious  damping  In  the  solution  (less  than  1  percent). 

Free  Field  Attenuation:  2-D  versus  3-D 

3.  The  displacement  versus  time  plot  for  the  three-dimensional  case 
shows  also  a  small  early  motion  starting  at  the  P-wave  arrival  time. 
Except  for  this  early  motion,  the  three  dimensional  displacement 
curve  retains  very  well  the  sinusoidal  shape  of  the  excitation 
because  both  positive  and  negative  peaks  a-e  almost  the  same  In 
amplitude,  and  the  displacement  after  a  loading  cycle  Is  perfectly 
zero  Instead  of  decreasing  asymtotlcally  as  in  the  two-dimensional 


case. 


4.  Damping  In  the  medium  not  only  reduces  the  amplitude  of  the  displace¬ 
ment  but  also  smoothens  out  all  abrupt  changes  In  the  motion  and 
results  In  a  very  smooth  displacement  curve.  The  velocity  of  the 
shear  wave  and  the  period  T  of  the  motion  appear  to  Increase  with 
damping.  These  effects  happen  In  both  the  two  and  the 
three-dimensional  cases. 

5.  Without  damping  In  the  medium,  and  at  distances  larger  than  twice  the 
wavelength,  the  log-log  plots  of  amplitude  versus  distance  are 
straight  lines  with  a  slope  equal  to  -1/2  and  -1  for  the  two  and  the 
three-dimensional  cases  respectively.  In  both  cases,  the  plots  are 
no  longer  a  straight  line  as  d/X  decreases.  The  maximum  curvature 
occurs  at  a  distance  of  approximately  one  wavelength  (d/X  =  1).  This 
deviation  from  a  straight  line  Is  due  to  the  combination  of  SV  and  P 
waves  caused  by  the  point  excitation. 

6.  In  a  medium  with  damping,  the  log-log  plots  of  amplitude  versus  dis¬ 
tance  are  not  straight  lines.  The  maximum  curvature  also  occurs 
around  d/X  *  1.  In  both  the  two  and  the  three-dimensional  cases  the 
curves  have  a  steeper  slope  as  damping  Increases  Indicating  that  the 
waves  attenuate  faster. 
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case. 
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Attenuation  in  a  Madium  with  an  Inclusion:  2-D  versus  3-D 

7.  For  both  the  two  and  the  three-dimensional  cases,  the  attenuation 
rate  with  the  presence  of  an  Inclusion  Is  virtually  the  same  as  that 
In  the  free  field  although  the  overall  amplitude  Is  reduced.  This 
Indicates  that  the  attenuation  rate  Is  not  a  good  Indicator  of  the 
presence  of  an  Inclusion. 

8.  The  reduction  In  amplitude  Is  more  pronounced  with  the  Increase  In 
the  size  of  an  Inclusion  (relative  to  the  wavelength).  In  the 
two-dimensional  case  without  damping,  the  reduction  Is  almost  con¬ 
stant  In  the  range  of  d/X  from  2  to  20.  As  the  damping  ratio  of  the 
medium  Increases,  the  effect  of  an  Inclusion  on  the  amplitude 
reduction  decreases  especially  at  larger  values  of  d/X.  For  all 
damping  ratios,  the  maximum  reduction  of  amplitude  occurs  at  a  dis¬ 
tance  of  about  1.25  times  the  wavelength  for  an  Inclusion  with  size 
equal  to  1/2  to  1/4  of  the  wavelength.  For  smaller  inclusions,  this 
distance  decreases  slightly.  In  the  three-dimensional  case,  the 
reduction  In  the  amplitude  due  to  an  Inclusion  Is  very  small  espe¬ 
cially  at  large  distances  and  for  large  damping  ratios,  which  Indi¬ 
cates  that  the  geometrical  damping  and  the  material  damping  play  a 
more  important  role  than  the  Inclusion  In  the  reduction  of  the  wave 
amplitude.  It  Is  thus  very  difficult  to  predict  the  existence  of  an 
Inclusion  from  the  amplitude  reduction  when  there  Is  Internal  damp¬ 
ing,  particularly  If  measurements  are  made  at  large  distances. 


ujjiEg?U4  'mwu'wmi 


176 


9.  Studies  conducted  for  the  two-dimensional  case  without  damping  show 
that,  for  different  stiffness  and  mass  density  of  the  Inclusion,  the 
rate  of  wave  attenuation  Is  practically  the  same.  An  Increase  In  the 
unit  mass  of  the  Inclusion  Increases  the  amplitude  slightly  at  all 
distance  ranges.  This  overall  Increase  In  amplitude  is  also 
achieved  with  a  decrease  In  the  stiffness  of  the  Inclusion. 

Effect  of  Inclusions  on  the  Wave  Velocity 

10.  The  apparent  velocity  of  the  propagating  waves  Is  affected,  as  could 
be  expected,  by  the  presence  of  an  Inclusion  along  the  wave  path.  It 
may  become  faster  or  slower  depending  upon  the  wave  velocity  (or 
stiffness)  of  the  Inclusion  relative  to  the  medium.  The  wave  veloci¬ 
ty  for  the  system  was  computed  from  the  first  arrival  time  of  the 
shear  wave  which  Is  determined  from  the  displacement  versus  time 
curve  at  the  point  where  the  major  wave  motion  starts.  This  point  Is 
easy  to  detect  at  some  distance  and  for  a  material  without  Internal 
damping  but  It  becomes  harder  to  locate  accurately  at  very  short  dis¬ 
tances  or  when  there  Is  material  damping. 

11.  The  maximum  effect  on  the  arrival  time  of  the  shear  waves  obviously 
occurs  when  the  travelling  path  of  the  wave  consists  largely  of 
Inclusions.  This  can  be  achieved  either  by  Increasing  the  size  or 
the  number  of  Inclusions,  or  more  specifically.  Increasing  the  total 
length  of  the  Inclusions  In  the  direction  of  the  wave  propagation. 


The  size  of  the  Inclusion  In  the  direction  perpendicular  to  that  of 
wave  propagation  also  affects  the  shear  wave  velocity. 

12.  The  solution  provided  by  one-dimensional  ray  theory,  computing 
directly  the  time  of  travel  In  both  media  gives  a  lower  bound  when 
the  Inclusion  Is  stlffer  than  the  surrounding  soil  and  an  upper  bound 
In  the  opposite  case. 

Effect  of  the  Inclusion  Shape 

Different  shapes  of  Inclusions  were  considered  In  this  study. 
They  were  classified  under  two  criteria:  same  length  of  the  principal 
wave  path;  and  same  total  area.  None  of  these  two  criteria  can  provide  a 
full  description  of  the  results  but  the  former  seems  to  be  the  better  of 
the  two. 

13.  Results  from  different  Inclusion  shapes  under  the  first  criterion, 
ranging  from  a  small  horizontal  rectangle  to  a  long  vertical  rectan¬ 
gle,  reveal  that  the  length  of  the  Inclusion  In  the  direction  perpen¬ 
dicular  to  the  principal  wave  path  also  affects  the  velocity  of  the 
wave.  The  results  show  that  the  longer  this  length  Is,  the  larger 
the  effect,  and  the  closer  the  curve  Is  to  the  ray  solution. 

14.  Results  for  the  inclusions  under  the  second  criterion,  for  a  hori¬ 
zontal  rectangle,  a  vertical  rectangle,  a  rhombus,  a  square,  and  a 
circle,  show  that  the  effect  on  the  the  arrival  time  of  the  wave  Is 


least  for  a  rhombic  shape  and  largest  for  a  square,  while  the  results 
from  the  circular  shape  fall  In  between.  The  horizontal  and  the  ver¬ 
tical  rectangles  also  show  results  that  are  close  to  those  for  a  cir¬ 
cular  shape.  These  two  cases,  however,  cannot  be  compared  well  with 
others  due  to  the  difference  In  the  length  of  the  wave  path  In  the 
Inclusions  (which  Is  used  as  the  normalizing  parameter). 

Multiple  Inclusion  Effect 

15.  Results  for  the  case  of  multiple  Inclusions  show  that  the  larger  the 
total  area  of  Inclusions  (or  the  number  of  Inclusions)  the  stronger 
the  effect  on  the  wave  velocity.  One  large  continuous  Inclusion 
affects  the  wave  velocity,  however,  more  than  a  cluster  of  several 
smaller  Inclusions  with  the  same  total  area. 

Effect  of  Inclusion  Stiffness 

16.  With  a  stlffer  inclusion  In  a  medium,  the  apparent  wave  velocity 
should  obviously  increase.  The  results  show  that  when  the  Inclusion 
is  less  stiff  than  the  medium,  the  change  in  the  arrival  time  Is  not 
as  pronounced  as  when  the  Inclusion  Is  stlffer  (than  the  medium). 
This  Is  because  the  wave  can  take  a  faster  route  around  the  Inclu- 
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Effect  of  tho  Unit  Mass  of  tho  Inclusion 

17.  For  the  square  Inclusion  considered  (a/X  =  1/8),  an  Increase  In  Its 
unit  mass  of  up  to  100  percent  does  not  significantly  change  the 
arrival  time.  The  Increase  In  the  arrival  time  Is  more  obvious  with 
a  larger  Increase  (above  300  percent)  In  the  unit  mass  of  the  Inclu¬ 
sion. 

This  work  was  Intended  to  demonstrate  the  numerical  formulation 
of  the  boundary  element  method  to  solve  the  problem  of  wave  propagation 
In  an  Infinite  medium  with  Inclusions  having  different  properties. 
Using  this  formulation,  further  parametric  studies  on  the  effect  of 
Inclusions  shoulo  be  conducted.  The  use  of  a  higher  order  element  to 
Improve  the  efficiency  of  the  model  should  also  be  investigated.  Another 
area  which  deserves  further  research  Is  that  of  wave  propagation  In  a 
nonhomogeneous  semi-infinite  medium.  This  may  be  done  either  by  discre¬ 
tizing  the  soil  surface  In  addition  to  the  soil-inclusion  contact  sur¬ 
face  (Dominguez,  1978),  or  by  using  the  dynamic  Green's  function  for  a 
semi-infinite  medium  (Johnson,  1974;  Apsel  and  Luco,  1983;  Luco  and 
Apsel,  1983).  In  this  way,  a  more  practical  model  for  wave  propagation 
In  soils,  Including  surface  waves  and  the  effects  of  Inclusions  can  be 
Investigated. 
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